# 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?

## 2 Answers

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.

Fortunately there are lots of literature on elements that work. Keywords are "Raviart-Thomas elements" and "Brezzi-Douglas-Marini elements". The monograph by Boffi-Brezzi-Fortin lists various finite element spaces for this particular problem.

In addition to @knl’s answer: you may want to check chapter 4 “Mixed Finite Elements” by Mardal (he is one of the developers of Dolfin which you mentioned) et al. of this book.

They discuss mixed Poisson problem, inf–sup stable FEs, and how it is implemented in Diffpack.