# The Laplacian

## 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{split} -\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{split}\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

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

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

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

 the keyword auto enables type inference, for more details see Wikipedia C++11 page.

Next 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.

auto Vh = Pch<2>( mesh );
auto u = Vh->element("u");
auto mu = doption(_name="mu");
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 g = expr( soption(_name="functions.g"), "g" );
auto v = Vh->element( g, "g" );
 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.

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);
a = integrate(_range=elements(mesh),
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");

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");

cout << "||u_h-g||_L2=" << normL2(_range=elements(mesh), _expr=idv(u)-g) << std::endl;

tic();
auto e = exporter( _mesh=mesh );
e->save();

toc("Exporter");
return 0;

}

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
//! 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);

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.

cout << "||u_h-g||_L2=" << normL2(_range=elements(mesh), _expr=idv(u)-g) << std::endl;

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
auto e = exporter( _mesh=mesh );
e->save();

## Testcases

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

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

### feelpp2d and feelpp3d

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

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

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