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));
%% Loop over quadrature points
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