I'm solving the Dirichlet problem for the Poisson equation in a 2d domain $D$:$$ \begin{cases}    \Delta u = 0 \quad \text{in $D$}, \\    u|_{\partial D} = u_0. \end{cases}$$I'm interested in finding the fluxes of the solution through the boundary and in order to get accurate fluxes (after a discussion in comments to this post) I decided to pass to a mixed formulation (see this fenics documentation page):$$ \begin{cases}   \nabla u = \sigma, \\   \nabla \cdot \sigma = 0, \\   u|_{\partial D} = u_0.  \end{cases}$$The weak form consists of two equations$$   \int_D (\sigma \cdot \tau + \nabla \cdot \tau \; u )dS = \int_{\partial D} \tau \cdot n \; u_0 \, ds \quad \forall \tau \in H^1(div), \\   \int_D \nabla \cdot \sigma \; v \; dS = 0 \quad \forall v \in L^2.$$Next, I consider a triangular finite element basis $\psi_i$ and I put $\sigma = (\sigma^i_x \psi_i, \sigma^j_y \psi_j)$, $\tau = (\tau^i_x \psi_i,\tau^j_y\psi_j)$, $u = u^i \psi_i$, $v = v^j \psi_j$ (where the summation signs are ommited). Putting these expressions in the above equations and taking into account that $\tau$, $v$ are arbitrary, I get the equations$$   \int_D \psi_i \psi_j dS \; \sigma_x^j + \int_D \partial_x \psi_i \psi_j dS \; u^j = \int_{\partial D} \psi_i n_x u_0 \, ds, \\   \int_D \psi_i \psi_j dS \; \sigma_y^j + \int_D \partial_y \psi_i \psi_j dS \; u^j = \int_{\partial D} \psi_i n_y u_0 \, ds, \\  \int_D \psi_i \partial_x \psi_j dS \; \sigma_x^j + \int_D \psi_i \partial_y \psi_j dS \; \sigma_y^j = 0,$$ where there is a summation over all indices $j$. Now I introduce the matrices$$  B_{ij} = \int_D \psi_i \psi_j \, dS, \; B_{x,ij} = \int_D \partial_x \psi_i \psi_j dS, \; B_{y,ij} = \int_D \partial_y \psi_i \psi_j \, dS,$$and the vectors$$   f_{x,i} = \int_{\partial D} \psi_i n_x u_0 \, ds, \; f_{y,i} = \int_{\partial D} \psi_i n_y u_0 \, ds.$$The discrete system can be formulated as$$  \left(\begin{array}{ccc}   0 & B_x^T & B_y^T \\   B_x & B & 0 \\   B_y & 0 & B \end{array}\right)\left(\begin{array}{c} u \\ \sigma_x \\ \sigma_y \end{array}\right)= \left(\begin{array}{c} 0 \\ f_x \\ f_y \end{array}\right).$$However, if I implement these formulas, I get some rubbish. Could you please point out what is wrong?
Mixed formulation of the Poisson equation (FEM)
finite element;pde;poisson