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?
$\begingroup$
$\endgroup$
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
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.
