[programming]

# Solving for the Laplacian

Feel++ allows you to solve a large range of partial differential equations. Here is one of the simplest

# The Laplacian

## 1. Problem statement

We are interested in this section in the conforming finite element approximation of the following problem:

Laplacian problem

Look for $u$ such that

\left\{\begin{aligned} -\Delta u &= f \text{ in } \Omega\\ u &= g \text{ on } \partial \Omega_D\\ \frac{\partial u}{\partial n} &=h \text{ on } \partial \Omega_N\\ \frac{\partial u}{\partial n} + u &=l \text{ on } \partial \Omega_R \end{aligned}\right.
 $\partial \Omega_D$, $\partial \Omega_N$ and $\partial \Omega_R$ can be empty sets. In the case $\partial \Omega_D =\partial \Omega_R = \emptyset$, then the solution is known up to a constant.
 In the implementation presented later, $\partial \Omega_D =\partial \Omega_N = \partial \Omega_R = \emptyset$, then we set Dirichlet boundary conditions all over the boundary. The problem then reads like a standard laplacian with inhomogeneous Dirichlet boundary conditions: Laplacian Problem with inhomogeneous Dirichlet conditions Look for $u$ such that Inhomogeneous Dirichlet Laplacian problem $-\Delta u = f\ \text{ in } \Omega,\quad u = g \text{ on } \partial \Omega$

## 2. Variational formulation

We assume that $f, h, l \in L^2(\Omega)$. The weak formulation of the problem then reads:

Laplacian problem variational formulation

Look for $u \in H^1_{g,\Gamma_D}(\Omega)$ such that

Variational formulation
$\displaystyle\int_\Omega \nabla u \cdot \nabla v +\int_{\Gamma_R} u v = \displaystyle \int_\Omega f\ v+ \int_{\Gamma_N} g\ v + \int_{\Gamma_R} l\ v,\quad \forall v \in H^1_{0,\Gamma_D}(\Omega)$

## 3. Conforming Approximation

We now turn to the finite element approximation using Lagrange finite element. We assume $\Omega$ to be a segment in 1D, a polygon in 2D or a polyhedron in 3D. We denote $V_\delta \subset H^1(\Omega)$ an approximation space such that $V_{g,\delta} \equiv P^k_{c,\delta}\cap H^1_{g,\Gamma_D}(\Omega)$.

Laplacian problem weak formulation

Look for $u_\delta \in V_\delta$ such that

$\displaystyle\int_{\Omega_\delta} \nabla u_{\delta} \cdot \nabla v_\delta +\int_{\Gamma_{R,\delta}} u_\delta\ v_\delta = \displaystyle \int_{\Omega_\delta} f\ v_\delta+ \int_{\Gamma_{N,\delta}} g\ v_\delta + \int_{\Gamma_{R,\delta}} l\ v_\delta,\quad \forall v_\delta \in V_{0,\delta}$
 from now on, we omit $\delta$ to lighten the notations. Be careful that it appears both the geometrical and approximation level.

## 4. Feel++ Implementation

In Feel++, $V_{g,\delta}$ is not built but rather $P^k_{c,\delta}$.

 The Dirichlet boundary conditions can be treated using different techniques and we use from now on the elimination technique.

We start with the mesh and the discretization setting by first defining `Vh=Pch<k>(mesh)` $\equiv P^k_{c,h}$, then elements of `Vh` and expressions `f`, `n` and `g` given by command line options or configuration file.

``````    tic();

tic();
auto Vh = Pch<2>( mesh );
auto u = Vh->element("u");
auto mu = expr(soption(_name="functions.mu")); // diffusion term
auto f = expr( soption(_name="functions.f"), "f" );
auto r_1 = expr( soption(_name="functions.a"), "a" ); // Robin left hand side expression
auto r_2 = expr( soption(_name="functions.b"), "b" ); // Robin right hand side expression
auto n = expr( soption(_name="functions.c"), "c" ); // Neumann expression
auto solution = expr( checker().solution(), "solution" );
auto g = checker().check()?solution:expr( soption(_name="functions.g"), "g" );
auto v = Vh->element( g, "g" );
toc("Vh");``````
 the keyword `auto` enables type inference, for more details see Wikipedia C++11 page.
 at the following line `` auto v = Vh->element( g, "g" );`` `v` is set to the expression `g`, which means more precisely that `v` is the interpolant of `g` in `Vh`.

the variational formulation is implemented below, we define the bilinear form `a` and linear form `l` and we set strongly the Dirichlet boundary conditions with the keyword `on` using elimination. If we don’t find `Dirichlet`, `Neumann` or `Robin` in the list of physical markers in the mesh data structure then we impose Dirichlet boundary conditions all over the boundary.

``````    tic();
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
_expr=f*id(v));
l+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_2*id(v));
l+=integrate(_range=markedfaces(mesh,"Neumann"), _expr=n*id(v));
toc("l");

tic();
auto a = form2( _trial=Vh, _test=Vh);
tic();
a = integrate(_range=elements(mesh),
toc("a");
a+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_1*idt(u)*id(v));
a+=on(_range=markedfaces(mesh,"Dirichlet"), _rhs=l, _element=u, _expr=g );
//! if no markers Robin Neumann or Dirichlet are present in the mesh then
//! impose Dirichlet boundary conditions over the entire boundary
if ( !mesh->hasAnyMarker({"Robin", "Neumann","Dirichlet"}) )
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );
toc("a");``````

We have the following correspondance:

Element sets Domain

`elements(mesh)`

$\Omega$

`boundaryfaces(mesh)`

$\partial \Omega$

`markedfaces(mesh,"Dirichlet")`

$\Gamma_D$

`markedfaces(mesh,"Neumann")`

$\Gamma_R$

`markedfaces(mesh,"Robin")`

$\Gamma_R$

next we solve the algebraic problem

Listing: solve algebraic system
``````    tic();
//! solve the linear system, find u s.t. a(u,v)=l(v) for all v
if ( !boption( "no-solve" ) )
a.solve(_rhs=l,_solution=u);
toc("a.solve");``````

next we compute the $L^2$ norm of $u_\delta-g$, it could serve as an $L^2$ error if $g$ was manufactured to be the exact solution of the Laplacian problem.

``````    // compute l2 and h1 norm of u-u_h where u=solution
auto norms = [=]( std::string const& solution ) ->std::map<std::string,double>
{
tic();
double l2 = normL2(_range=elements(mesh), _expr=idv(u)-expr(solution) );
toc("L2 error norm");
tic();
toc("H1 error norm");
return { { "L2", l2 }, {  "H1", h1 } };
};

int status = checker().runOnce( norms, rate::hp( mesh->hMax(), Vh->fe()->order() ) );``````

and finally we export the results, by default it is in the ensight gold format and the files can be read with Paraview and Ensight. We save both $u$ and $g$.

Listing: export Laplacian results
``````    tic();
auto e = exporter( _mesh=mesh );
if ( checker().check() )
{
v.on(_range=elements(mesh), _expr=solution );
}
e->save();
toc("Exporter");``````

## 5. Testcases

The Feel++ Implementation comes with testcases in 2D and 3D.

### 5.1. circle

`circle` is a 2D testcase where $\Omega$ is a disk whose boundary has been split such that $\partial \Omega=\partial \Omega_D \cup \partial \Omega_N \cup \partial \Omega_R$.

Here are some results we can observe after use the following command

``````cd Testcases/quickstart/circle
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file circle.cfg``````

This give us some data such as solution of our problem or the mesh used in the application.

 Solution $u_\delta$ Mesh

### 5.2. feelpp2d and feelpp3d

This testcase solves the Laplacian problem in $\Omega$ an quadrangle or hexadra containing the letters of Feel++

#### 5.2.1. feelpp2d

After running the following command

``````cd Testcases/quickstart/feelpp2d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file feelpp2d.cfg``````

we obtain the result $u_\delta$ and also the mesh

 /images/Laplacian/TestCases/Feelpp2d/meshfeelpp2d.png[] Solution $u_\delta$ Mesh

#### 5.2.2. feelpp3d

We can launch this application with the current line

``````cd Testcases/quickstart/feelpp3d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_3d --config-file feelpp3d.cfg``````

When it’s finish, we can extract some informations

 Solution $u_\delta$ Mesh