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

``    auto mesh = loadMesh(_mesh=new Mesh<Simplex<FEELPP_DIM,1>>);``
 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