# Solving Poisson equations as mixed Laplace using $RT_0-P_0$ pair

I'm trying to solve $$\begin{cases} - \Delta p=f \text{ in } \Omega\\ p=0 \text{ on } \partial \Omega \end{cases}$$ with in $$\Omega = [-1,1]^2$$ by writing it as

$$\begin{cases} u + \nabla p=0 \\ -\operatorname{div}(u) = -f \\ p = 0 \text{ on } \partial \Omega \end{cases}$$

whose weak form is $$\begin{cases}(v,u) - (\operatorname{div}(v),p) = 0 \qquad \forall v \in V\\ -(\operatorname{div}(u),q) = -(f,q) \qquad \forall q \in Q\end{cases}$$

where $$V=H^{\operatorname{div}}(\Omega)$$ and $$Q=L^2(\Omega)$$. To solve it, I decided to use the inf-sup stable couple $$V_h=RT_0$$ (for the velocity) and $$Q=P_0$$ for the pressure. The basis functions for $$RT_0$$ in the reference triangle $$\hat{K}$$ are $$\hat{\phi_1} = \sqrt{2}(\hat{x},\hat{y})$$ $$\hat{\phi_2} = (-1+\hat{x},\hat{y})$$ $$\hat{\phi_3} = (\hat{x},-1+\hat{y})$$

In terms of finite element matrices, we have a saddle point problem and the element matrices $$A^K$$ and $$B^K$$ ($$K$$ is a triangle) have components:

$$a_{ij}^K = \frac{1}{|\det(B_K)|}\int_{\hat{K}} [\text{sign}_i^K][ \text{sign}_j^K] B_K \hat{\phi_i} \cdot B_K \hat{\phi_j}$$

$$b_j^K=-\frac{1}{|\det(B_K)|} \int_{\hat{K}} [\text{sign}_j^K] \operatorname{div}(\hat{\phi_j})$$

where $$B_K$$ is the matrix in the classical affine mapping $$F_K:\hat{K} \rightarrow K$$, $$F_K(\hat{\boldsymbol{x}}) = B_K \hat{\boldsymbol{x}} + \boldsymbol{b_K}$$.

I've been implementing this in MatLab for two days, but the condition number of the whole saddle point system is infinite.

• The boundary condition $$p=0 \text{ on } \partial \Omega$$ should be weakly imposed, so I did not change the matrix after the assemble() function. If that is correct, then the problem must be inside my assemble() function, in particular in the distribution of the entries of $$B$$. Since I have $$1$$ DoF per triangle for the pressure, I have a 3x1 vector for each element $$K$$.

I think the following code is really "didactic":

• here the inputs p,t are the result of the MatLab function initmesh, and force is a function handle with the forcing term.

• RT_shapes is a function that evaluates at the points the RT basis functions and also computes the divergence for each function (which happens to be a constant vector $$[2 \sqrt{2},2,2]$$)

• In my tests, I am assuming $$f$$ s.t. the solution is $$(x^2-1)(y^2-1)$$, which indeed satisfies homogeneous Dirichlet.

Do you spot any error in my reasoning? I really don't see it, and any hint is really welcome!

function [A,B,F] = assemble(p,t,force)
[rspoints,qwgts] = GaussPoints(4);
np = size(p,2); %N points
nt = size(t,2); %N elements
A = sparse(np,np); %N_DoFs x N_DoFs
B = sparse(nt,np); % N_triangles x N_DoFs
F = zeros(nt,1);

for K=1:nt
l2g = t(1:3,K); %global node indices for element K
tmp = l2g([2 3 1]) - l2g([3 1 2]);
signs = tmp ./ abs(tmp);

x = p(1,l2g); %x coords
y = p(2,l2g); %y coords

BK = [x(2)-x(1), x(3)-x(1);y(2)-y(1),y(3)-y(1)];
bK = [x(1);y(1)];
detBK = det(BK);

detBK_inv = 1/(abs(detBK));
for q=1:length(qwgts)
r = rspoints(q,1); %x coordinate q-th quadrature point
s = rspoints(q,2); %y coordinate q-th quadrature point
[phi,divphi] = RT_shapes(r,s);

JxW=qwgts(q)*detBK/2.0;
physical_coords = BK*[r;s] + bK; %physical coordinates of the current quadrature points
xp = physical_coords(1);
yp = physical_coords(2);

val_rhs = -force(xp,yp)*1.0*JxW;
F(K) = F(K) + val_rhs;

for i=1:3
for j=1:3
val_A =  signs(i)*signs(j)*detBK_inv* dot(BK*phi(:,i),BK*phi(:,j))*qwgts(q);
A(l2g(i),l2g(j)) = A(l2g(i),l2g(j)) + val_A;

end
val_B = - signs(i)* detBK_inv*divphi(i)*qwgts(q);
B(K,l2g(i)) = B(K, l2g(i)) + val_B;

end

end

end


# EDIT

@knl spotted a fatal typo in my code above, i.e. the indexing given by l2g was the one referred to vertex, not to triangle edges. I used a suitable routine, found in the appendix of the book by Larson-Bengzon, named Tri2Edge that numbers the edges of a triangle mesh.

Now, the following code is the last version.

• I'm using the indices of the edges to decide which entry of the matrices $$A$$ and $$B$$ has to be filled.

• Let $$N_t$$ the number of elements, and $$N_e$$ the number of edges. The matrix $$A$$ is a $$N_e \times N_e$$, while $$B$$ is a $$N_e \times N_t$$. The $$F$$ term in the rhs has size $$N_t$$.

• Again, homogeneous Dirichlet BC are assumed.

• Crucially, I noticed that in my reference paper for the implementation, the formula $$(8)$$ for $$b_j^K$$ is not what they implemented, since if you go to the first snippet, you may see that they multiplied by detJ the entries of $$B$$, while in that formula they divide their divergence by detJ. I don't know why, but if I do not multiply by detJ, I obtain the correct solution, as can be seen by the following graph for the pressure, and the correct $$L^2$$ convergence for the pressure (order $$1$$, as I am using 1 DoF per triangle).

[![enter image description here][2]][2]

function [A,B,F] = assemble(p,t,t2e,force)
[rspoints,qwgts] = GaussPoints(4);
nt = size(t,2); %N_triangles
ne = max(t2e(:)); %N_edges
A = sparse(ne,ne); %N_edges x N_edges
B = sparse(nt,ne); % N_triangles x N_edges
F = zeros(nt,1); %N_triangles

for K=1:nt
l2g = t(1:3,K); %global node indices for element K
edges = t2e(K,:);%global edges indices for element K
tmp = l2g([2 3 1]) - l2g([3 1 2]);
signs = tmp ./ abs(tmp);

x = p(1,l2g); %x coords
y = p(2,l2g); %y coords

BK = [x(2)-x(1), x(3)-x(1);y(2)-y(1),y(3)-y(1)];
bK = [x(1);y(1)];
detBK = det(BK);

detBK_inv = 1/abs(detBK);
for q=1:length(qwgts)
r = rspoints(q,1); %x coordinate q-th quadrature point
s = rspoints(q,2); %y coordinate q-th quadrature point
[phi,divphi] = RT_shapes(r,s);

JxW=qwgts(q)*detBK;
physical_coords = BK*[r;s] + bK; %physical coordinates of the current quadrature points
xp = physical_coords(1);
yp = physical_coords(2);

val_rhs = -force(xp,yp)*1.0*JxW;
F(K) = F(K) + val_rhs;

for i=1:3
for j=1:3
val_A =  detBK_inv*dot(signs(i)*BK*phi(:,i),signs(j)*BK*phi(:,j))*qwgts(q);
A(edges(i),edges(j)) = A(edges(i),edges(j)) + val_A;

end
val_B = -signs(i)*divphi(i)*qwgts(q);
B(K,edges(i)) = B(K, edges(i)) + val_B;

end

end

end

• What is the nullity of your global matrix? How do you make sure that the BC for the pressure is satisfied (in the weak sense)? Dec 31, 2021 at 19:38
• Before we run code for RT/Hdiv space, we run a unit test, checking if the volume integral divergence is equal to the flux on the boundary (.i.e. Gauss's theorem). That will validate your base functions/integration and Piola transform. We set the coefficient of every base to one. Our code is generic for arbitrary base order, see mofem.eng.gla.ac.uk/mofem/html/… Dec 31, 2021 at 20:49
• I found a MATLAB implementation: mathematik.hu-berlin.de/~cc/download/public/software/… Dec 31, 2021 at 22:22
• @ZoltánCsáti The nullity for $B'$ returns a matrix with 3 columns. I know that it should return a single vector $1$, or a scaled copy of $1$, so the problem must be there. For what concerns the pressure, I was assuming that the condition $p=0$ on $\partial \Omega$ is automatically satisfied since the term on the weak form has been eliminated Jan 1 at 11:06
• @ZoltánCsáti I've checked everything again and I can't spot the problem. My snippet wants just to compute the local contributions to the mass matrix $A$ and the "divergence" $B$. The resulting global matrix should not be post-processed, since I want homogeneous Dirichlet on the whole domain. I'm really puzzled honestly Jan 2 at 17:33

You need to use a different kind of indexing in l2g. Right now you are using the indices of the global vertices when deciding which rows and columns to fill. This causes wrong entries to be summed when building the global matrix. You need to instead uniquely index all edges of your mesh and use those indices in l2g.

• Honestly, I just found the expression for those functions and used them. I took them from siam.org/Portals/0/Publications/SIURO/Vol12/S01743.pdf (page 250) Jan 1 at 10:59
• Alright, in the DOF you are dividing by the length of the edge. Then it's fine I guess.
– knl
Jan 1 at 15:27
• RT0 basis functions are one per edge, not one per vertex. The indexing is different.
– knl
Jan 2 at 22:05
• Yes, we are applying two transformations. $F_K$ defines the shape of the global triangle. Piola makes sure that the derivative DOFs are preserved through the transformation. You can think about it in 1D. If you have $f(u) = u(a)$ as a DOF we can transform basis function as $\psi(x) = \hat{\psi}(\hat{x})$ and, as a consequence, $f(\hat{\psi}) = f(\psi)$. However, if you have $g(u) = u^\prime(a)$ as a DOF we would have $g(\hat{\psi}) \not= g(\psi)$. Hence, we cannot do $\psi(x) = \hat{\psi}(\hat{x})$ and need to invent something else. (Piola preserves normal derivatives in 2D and 3D).
– knl
Jan 5 at 7:31
• Note that even in standard $H^1$ FEM we are doing two transformations, the other one is just so trivial that you never think about it as a transformation. ($\psi(x) = \hat{\psi}(\hat{x})$)
– knl
Jan 5 at 7:34