# Finite elements convergence issue with 2D elliptic equation

I deal with a system of coupled 2D Helmholtz-like equations solved via the P1 FEM on a given geometry. Let's consider, for instance, the following simplified coupled problem:

for $$i\in[-I,I]$$ we have $$\Delta a_i+k^2\sum_{j\in[-I_1,I_1]}c_ja_j=0$$ where $$[-I_1,I_1]\subset[-I,I]$$ and $$c_i\in\mathbb{C}$$. Only one of these equation has non zero Dirichlet condition at some node $$a_i(n_s)=1$$ for $$i=0$$ (source). The others ($$i\neq0$$) are either defined with homogeneous Dirichlet, or ABC (transparent).

The final solution is then comprised of the family $$\{a_i\}_i$$. My issue is the following. We have as expected $$a_0(n_s)=1$$, and the global shape of $$a_0$$ is more or less as usual (concentric circles when full ABC), but its amplitude decreases when we refine the mesh (while keeping, as expected, $$a_0(n_s)=1$$), without apparent convergence. For the other solutions, the shape is constant and exactly as expected, but the amplitude decreases linearly when we refine the mesh. This behavior is very explicit, and I expect that it comes from a well known rough mistake, but up to know I have not found it.

As I think this issue rises from the way to code the assembly, I think it also appears for one single 2D Helmholtz solver, as toy model, and for which I edit the following code below. Not very fast (the loops are normally parallelized), I wrote it explicit to make it easy to debug (at least I hope).

Thanks in advance for any help!!

Create the mesh

model   = createpde;
Domain  = [3 4 0 1 1 0 0 0 1 1]';  % Unite square domain
g       = decsg(Domain,'D',char('D')');
geometryFromEdges(model,g);
mesh    = generateMesh(model,'Hmax',0.01,'Hmin',0.001,'GeometricOrder','Linear');
[p,e,t] = meshToPet(model.Mesh);   % Triangulation

t = t(1 : 3, :).';   % Get column PET format with usefull labels
p = p.';
e = e([1 2 5],:)';

Antenna = [0.5 0.5]; % Expected position of the source
T = triangulation(t(:,1 : 3),p); % Triangulation data from mesh


Indices of physical regions. We define which boundary supports ABC or homogeneous Dirichlet condition. Set the antenna position.

id_E1 = find(e(:,3) == 1);  % Get edge indices of each of the 4 sides of the domain
id_E2 = find(e(:,3) == 2);
id_E3 = find(e(:,3) == 3);
id_E4 = find(e(:,3) == 4);

id_ABC    = [id_E1; id_E2; id_E3; id_E4];  % Choose edges supporting Absorbing Boundary Condition
Dir_guess = [];                % Choose edges supporting homogeneous Dirichlet condition

id_Dir      = unique(e(Dir_guess,1:2));   % Extract indices of the nodes supporting homogeneous Dirichlet boundary conditions (reflecting parts)
id_Antenna  = nearestNeighbor(T,Antenna);               % Index of the node where we apply inhomogeneous Dirichlet: i.e. the antenna (found via nearestNeighbor)
id_DirNodes = [id_Antenna; setdiff(id_Dir,id_Antenna)]; % Indices of all the Dirichlet nodes (Antenna + reflecting parts)

nb_t    = size(t,1);               % Number of elements
nb_p    = size(p,1);               % Number of nodes
nb_eABC = length(id_ABC);          % Number of ABC edges
nb_DirNodes = length(id_DirNodes); % Number of Dirichlet nodes

id_DoF = setdiff((1 : nb_p).',id_DirNodes);  % Indices of the free nodes (ABC nodes + interior nodes)


Physical parameters: wavenumber and amplitude of emission.

k = 30;                % Wavenumber of the Helmholtz equation
Source_Amplitude = 1;  % Amplitude of emission


FEM assembly. Set tables $$tM$$ (for elementary mass), $$tK$$ (elementary stiffness), and $$tbM$$(elementary edge mass) with the elementary matrices at each column: minimal cost for storage. For the element $$T_r$$ of the covering of the domain, given the basis $$\lbrace \phi_j\rbrace_{j\leq nb_p}$$ ($$nb_p$$ number of nodes) of linear shape functions being 1 on node $$j$$. Then the elementary matrices write $$M^{i,j}_r=k^2\int_{T_r}\phi_i.\phi_j dT_r,$$ $$K^{i,j}_r=-\int_{T_r}\nabla\phi_i.\nabla\phi_j dT_r$$ with $$1\leq i,j\leq3$$. Given the absorbing boundary condition, the IPP yields for edge $$E_l$$ $$bM^{i,j}_l=ik\int_{E_l}\phi_i.\phi_j dE_l,$$ with $$1\leq i,j\leq2$$. And basically I evaluate them using the areas and lengths of the involved elements, directly as shown in the code. For me the only Jacobian to evaluate is the one from the reference unit element to the one we consider, which requires the determinant.

tM  = zeros(9,nb_t);   % Initialize tables of elementary matrices (M: mass, K: stiffness, i-th column = reshaped elementary matrix of triangle i)
tK  = zeros(9,nb_t);
tbM = zeros(9,nb_t);   % Initialize table of edge mass elementary matrices in (explain)

for i = 1 : nb_t  % Loop of the triangles
%---Characterization of triangular element i
nodes_t = t(i,1:3);                     % Indices of the 3 nodes of triangle i
Pe      = [ones(3,1), p(nodes_t,1:2)];  % 3x3 matrix with lines (1 Xcorner Ycorner)
CoefT   = inv(Pe);                      % Columns of CoefT are coefs a,b,c in shape function phi = a+bx+cy
Area    = abs(det(Pe))/2;
%---Elementary stiffness of element i
grad = CoefT(2:3,:);            % P1: grad = the coefs
Ke   = Area * (grad.' * grad);  % Elementary stiffness
Ke   = reshape(Ke,9,1);         % Reshape in 9x1 to load in tK
%---Elementary mass of element i
GaussPt = [1/2, 1/2, 0; 0, 1/2, 1/2; 1/2, 0, 1/2] * Pe(:,2 : 3); % Location of gauss points on ref triangle
MeanPt  = [ones(3,1) GaussPt] * CoefT;                  % Apply to shape function
Me      = Area/3 * (MeanPt(1,:).' * MeanPt(1,:) + ...
MeanPt(2,:).' * MeanPt(2,:) + MeanPt(3,:).' * MeanPt(3,:)); % Mass elementary matrix: jac * (Gauss quad of phi * phi)
Me = reshape(Me,9,1);    % Reshape in 9x1 to load in tM
%---Load the elementary matrices on the tables. It seems that this is a wrong way to proceed
tM(:,i) = Me;
tK(:,i) = Ke;
end


For simplicity, the table $$tbM$$ is rewritten in the format of the other tables, using the connectivity lists.

I  = t(:,[1 2 3 1 2 3 1 2 3]).';  % Connectivity list of triangles, rows
J  = t(:,[1 1 1 2 2 2 3 3 3]).';  % Connectivity list of triangles, columns

Ib = e(id_ABC(:),[1 2 1 2]).';    % Same for lines (edge elements)
Jb = e(id_ABC(:),[1 1 2 2]).';

ID = edgeAttachments(T,e(id_ABC,1),e(id_ABC,2)); % From each considered edges, find associated triangle

for i = 1 : nb_eABC % Loop over the ABC edges
%---Elementary edge mass matrix bMe of edge i
nodes_e = e(id_ABC(i),1 : 2);  % Indices of the 2 nodes of edge i
Le      = norm(p(nodes_e(2),1 : 2) - p(nodes_e(1),1 : 2)); % Length of edge i
bMe     = Le/6 * [2; 1; 1; 2]; % Resolution of 1D integral
%---Spread contributions of bMe, according to triangle format of the previous tables
idTRIelem = [I(:,ID{i}) J(:,ID{i})]; % Find connectivity list of triangle linked to edge i
iE1  = [Ib(1,i) Jb(1,i)];            % Connectivity of each contribution of the elementary matrices
iE2  = [Ib(2,i) Jb(2,i)];
iE3  = [Ib(3,i) Jb(3,i)];
iE4  = [Ib(4,i) Jb(4,i)];
Pos1 = find(idTRIelem(:,1) == iE1(1) & idTRIelem(:,2) == iE1(2)); % Find the corresponding position in the global table
Pos2 = find(idTRIelem(:,1) == iE2(1) & idTRIelem(:,2) == iE2(2));
Pos3 = find(idTRIelem(:,1) == iE3(1) & idTRIelem(:,2) == iE3(2));
Pos4 = find(idTRIelem(:,1) == iE4(1) & idTRIelem(:,2) == iE4(2));
%---Load the contributions in tbM
tbM(Pos1,ID{i}) = tbM(Pos1,ID{i}) + bMe(1);
tbM(Pos2,ID{i}) = tbM(Pos2,ID{i}) + bMe(2);
tbM(Pos3,ID{i}) = tbM(Pos3,ID{i}) + bMe(3);
tbM(Pos4,ID{i}) = tbM(Pos4,ID{i}) + bMe(4);
end


Build the global matrix: global table $$tA$$ and sparsification under matrix format.

tA    = k^2 * tM + 1i * k * tbM - tK;         % Table with elementary matrices of global system
A     = sparse(I,J,tA,nb_p,nb_p); % Assemble sparse global FEM matrix
A_DoF = A(id_DoF,id_DoF);  % Extract the DoF for the resolution (Dirichlet nodes sent to right hand side)


Right hand side with Dirichlet conditions

DirBound = [Source_Amplitude; zeros(nb_DirNodes - length(Source_Amplitude),1)];    % Set Dirichlet condition vector
b        = -A(id_DoF,id_DirNodes) * DirBound;     % Right hand side, enforcing the Dirichlet conditions


Resolution of the linear system

U_vect = A_DoF\b;


Complete the solution with boundary nodes

U              = zeros(nb_p,1);
U(id_DirNodes) = DirBound;  % Fill with nodes with Dirichlet condition (Antenna + reflecting parts)
U(id_DoF)      = U_vect;    % Fill with DoF (interior nodes + ABS edge nodes)


Plot

trisurf(t(:,1:3), p(:,1), p(:,2),real(U),'facecolor','interp');
set(0,'DefaultFigureColormap',jet()); shading interp;

• There's clearly too much code here for anyone to say for sure. Have you tried with a simple problem where you have only two of these functions and where you may be able to write down the exact solution? – Wolfgang Bangerth Dec 29 '20 at 5:54
• You mean two solutions $a_i$? Basically yes, and this code is even to solve the more simple case of only one solution, from the single helmholtz equation – DgSl Dec 29 '20 at 10:13
• So then verify that your code is correct for one equation. If that's correct, go to the two-equation case, etc. – Wolfgang Bangerth Dec 30 '20 at 23:12
• I verified, also using the exact solution in 2d with hankel's functions, it converges, but quite slowly. I guess that P1 elements are not accurate, and then converge slowlier that quadratic elements, right? – DgSl Dec 31 '20 at 22:56
• Yes, linear elements are slower to converge. But you need to say what you mean by convergence: Convergence as a function of mesh size? Convergence as a function of nonlinear iterations (if you have any)? Etc. – Wolfgang Bangerth Jan 1 at 23:43