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

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).

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);
    %% 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;
        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