# Mixed formulation of the Poisson equation (FEM)

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?

For this problem you cannot use arbitrary finite elements. The finite element spaces should satisfy the so-called discrete inf-sup condition. You did not mention what elements you used so I assume that you used something that does not work like $P_1$ for each variable.