The Poisson equation with Coefficient Form PDEs toolbox

We consider a conformal approximation of problem by Lagrange finite elements. We assume for simplicity that \(\Omega\) is a \(\mathbb{R}^2\) polygon or a \(\mathbb{R}^3\) polyhedron. We consider a regular and conformal family of affine meshes of \(\Omega\) which we denote \(\left\{\mathcal{T}_h\right\}_{h>0}\).

We choose as finite element of reference \(\left\{\widehat{K}, \widehat{P}, \widehat{\Sigma}\right\}\) a finite element of Lagrange such that \(\mathbb{P}_k \subset \widehat{P}\) and \(k+1>\frac{d}{2}\)

We pose

where \(T_K\) is the geometric transformation sending \(\widehat{K}\) into \(K\). The space \(P_{\mathrm{c}, h}^k\) is the space of the functions of \(\mathcal{C}^0(\bar{\Omega})\) which are continuous on the boundary of \(\Omega\) and which are piecewise polynomials of degree \(k\) on each element of \(\mathcal{T}_h\).

In order to construct an approximation space which is included in \(V=H_0^1(\Omega)\), we pose

The elements of \(V_h\) are the functions of \(P_{\mathrm{c}, h}^k\) which are \(0\) on the boundary of \(\Omega\).

We consider the following approximate problem:

which is clearly well-posed since \(a\) is coercive on \(V\) and \(V_h \subset V\).

We now consider an example to illustrate the previous theorem.

We start by initializing the Feel++ environment.

Then we consider the \(\mathbb{R}^2\) domain \(\Omega\) defined by \(\Omega=\left\{x \in \mathbb{R}^2; 0 \leq x_1 \leq 1, 0 \leq x_2 \leq 1\right\}\). We consider the following mesh \(\mathcal{T}_h\) of \(\Omega\) with \(h=0.1\).

Then we consider the following right hand side \(f\) and the exact solution \(u(x,y) = \sin(2 \pi x) \sin(2 \pi y)\) such that

that way we can check the Convergence theorem of the conforming approximation with the exact solution.

The Coefficient Form PDEs toolbox allows to solve the problem problem with the following code:

We can proceed with the visualisation of the field u using the following code using pyvista:

In 3D, the plots are more difficult to visualise. We can use the pyvista package to visualise the mesh and the field u:

Then we can run the following code to:

In 3D, we use the same code to verify the convergence rates of the error in \(L^2\) and \(H^1\) norms. However we do this for different mesh sizes and only in \(\mathbb{P}^1\) to reduce the computational cost which increases significantly with the dimension.

Again, we observe that the proper convergence rates are obtained.

Finally we postprocess the results using plotly. The following code allows to

and now we plot the convergence rates in \(L^2\) and \(H^1\) norms for the 3D case.

Given a function \(f \in L^2(\Omega)\) and a function \(g \in \mathcal{C}^{0,1}(\partial \Omega)\) ( \(g\) is lipschitzian on \(\partial \Omega\)), we look for a function \(u: \Omega \rightarrow \mathbb{R}\) such that

The hypothesis \(g \in \mathcal{C}^{0,1}(\partial \Omega)\) allows us to assert that there exists a lifting \(u_g\) of \(g\) in \(H^1(\Omega)\), that is to say that there exists a function \(u_g\) in \(H^1(\Omega)\) such that \(\left.u_g\right|_{\partial \Omega}=g\).

Under these conditions, we perform the change of unknown \(u_0=u-u_g\) and consider the following weak formulation:

By the Lax-Milgram lemma, this problem is well posed. We are interested in a conformal approximation of the problem by Lagrange finite elements. We use the discrete framework described previously. We suppose that the data \(g\) is regular enough to admit a lifting \(u_g\) in \(\mathcal{C}^0(\bar{\Omega}) \cap H^1(\Omega)\). We denote \(\mathcal{I}_h^{\mathrm{Lag}}\) the interpolation operator associated with the mesh \(\mathcal{T}_h\) and the finite Lagrangian element of reference \(\widehat{K}, \widehat{P}, \widehat{\Sigma}}\). Recall that \(P_{\mathrm{c}, h}^k\) denotes the \(H^1\)-conformal space based on this finite element and that \(V_h\) is the \(H_0^1\)-conformal approximation space defined in (5.9).

Let \(N=\operatorname{dim} P_{\mathrm{c}, h}^k\) be given. We denote by \(\left\{\varphi_1, \ldots, \varphi_N\right\}\) the nodal basis of \(P_{\mathrm{c}, h}^k\) and by \(\left\{a_1, \ldots, a_N\right\}\) the associated nodes. By definition, for \(u \in \mathcal{C}^0(\bar{\Omega})\), we have

and we also introduce the surface Lagrange interpolation

We consider the approximated problem

which is clearly well posed. We pose \(u_h=u_{0 h}+\mathcal{I}_h^{\mathrm{Lag}} u_g\) so that \(\left.u_h\right|_{partial \Omega}\) coincides with the Lagrange interpolated surface of \(g\).

\(u_h\) is the solution of the problem

Given a strictly positive real \(\lambda\), a function \(f \in L^2(\Omega)\) and a function \(g \in L^2(\partial \Omega)\), we look for a function \(u: \Omega \rightarrow \mathbb{R}\) such that

where \(\partial_n u\) denotes the normal derivative of \(u\) on the boundar.

The weak formulation of the Neumann problem reads

where

and

The well-posedness of the weak formulation results from the coercivity of the bilinear form \(a_{\mathrm{N}} \operatorname{on} H^1(\Omega)\). The problem can be approximated by Lagrange finite elements. We take the discrete framework described for the Conforming Approximation. We consider the following approximated problem:

which is clearly well-posed.

Finally, the convergence analysis of the approximatio problem leads, under the assumptions of Theorem, to the same estimates as for the homogeneous Dirichlet problem.

We can also consider, at the cost of some technical difficulties, the Neumann problem with \(\lambda=0\). We look for a function \(u: \Omega \rightarrow \mathbb{R}\) such that

A solution of the problem with the condition is only determined to within one additive constant. We should then look for the solution in the functional space

The weak formulation is the following:

with the bilinear form $a$ such that \(a(v, w)=\int_{\Omega} \nabla v \cdot \nabla w\).

The well-posedness of problem results from the coercivity of the bilinear form a on \(H_*^1(\Omega)\), this last property being itself a consequence of the inequality below.

The problem can be approximated by Lagrangian finite elements, which leads to the following approximated problem:

By noting this linear system in the form

and introducing the vector \(Z=(1, \ldots, 1)^T \in \mathbb{R}^N\) where \(N=\operatorname{dim} P_{\mathrm{c}, h}^k\), we find that \(\operatorname{Ker}\left(\mathcal{A}^*\right)=\operatorname{span}(Z), \operatorname{Im}\left(\mathcal{A}^*\right)=Z^{\perp}\) and \(F^* \in Z^{\perp}\). Therefore, the linear system has infinitely many solutions. One of them can be approximated using an iterative method, for example the conjugate gradient method. We denote by \(U^{\circ}\) the approximation provided by the iterative method. The solution of the linear system whose components have zero mean is obtained by posing

where \((\cdot, \cdot)_{\mathbb{R}^N}\) denotes the Euclidean scalar product on \(\mathbb{R}^N\). \ The postprocessing of the solution removes the mean value of the solution \(U^{\circ}\) and is necessary to obtain the zero mean solution of problem.

An alternative way to obtain the solution of linear system is to use the method of Lagrange multipliers to impose the condition \(\int_{\Omega} u=0\).

Given a function \(f \in L^2(\Omega)\) and a function \(g \in L^2(\partial \Omega)\), we look for a function \(u: \Omega \rightarrow \mathbb{R}\) such that

where \(\gamma_R>0\) is a real parameter.

The weak formulation of the problem reads

where the linear form \(f_{\mathrm{N}}\) is defined here and the bilinear form \(a_{\mathrm{R}}\) is such that

The well-posedness of weak problem results from the coercivity of \(a_{mathrm{R}}[\) on \(H^1(\Omega)\), this last property being itself a consequence of the fact that there exists a constant \(\rho_{\Omega}\) such that

Starting from the initial discrete framework, we are interested in an approximation by Lagrange finite elements. We consider the following approximated problem:

which is clearly well-posed.

Finally, the convergence analysis for the \(u_h\) leads, under the assumptions of <<thm:1,the theorem>, to the same estimates as for the homogeneous Dirichlet problem.

Consider a differential operator of the form ]

where $\sigma, \beta$ and $\mu$ are functions defined on $\Omega$ and with values in $\mathbb{R}^{d, d}, \mathbb{R}^d$ and $\mathbb{R}$ respectively. The operator is used, for example, in the modeling of advection-diffusion-reaction problems; the first term of the right-hand side is the diffusive term, the second the advective term and the third the reactive term.

In the following, we assume that \(\sigma \in\left[L^{\infty}(\Omega)\right]^{d, d}, \beta \in\left[L^{\infty}(\Omega)\right]^d, \nabla \cdot \beta \in L^{\infty}(\Omega)\) and \(\mu \in L^{\infty}(\Omega)\). Moreover, we assume that the operator \(\mathcal{L}\) is elliptic in the following sense.

Given a function \(f \in L^2(\Omega)\), we look for a function \(u: \Omega \rightarrow \mathbb{R}\) such that

The weak formulation of the strong formulation is the following:

with the bilinear form \(a_{\sigma \beta \mu}\) such that for all \((v, w) \in H_0^1(\Omega) \times H_0^1(\Omega)\),

We pose \(p=\inf {\operatorname{ess}}{ }_{x \in \Omega}\left(\mu-\frac{1}{2} \nabla \cdot \beta\right)\) and we recall that \(\ell_{\Omega}\) is the constant intervening in the Poincaré inequality above.

Then, we show that under the condition

the bilinear form \(a_{\sigma \beta \mu}\) is coercive on \(H_0^1(\Omega)\), so that, thanks to the Lax-Milgram lemma, the problem is well posed. The condition being a minimization of \(\sigma_0\), we retain that the coercivity is guaranteed provided that \(\sigma_0\) is sufficiently large, i.e. in the dominant diffusion regime.

The conformal approximation by Lagrangian finite elements uses the discrete framework introduced for the Laplacian. Finally, one can consider instead of non-homogeneous Dirichlet, Neumann or Robin boundary conditions.