:leveloffset:+1

# 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

# Levelset

Having the possibility to determine where two regions meeting can be really useful in some scientific domains. That’s why the levelset method is born.

## Levelset introduction

### Levelset function

By using a scalar function \phi, define on all regions as the null value is obtained when it’s placed on an interface of two domains.

We denote \Omega_1 and \Omega_2 two domains with \Gamma the interface betwen them. Then \phi can be define as

\phi(\boldsymbol{x}) = \left\{ \begin{array}{cccc} \text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in &\Omega_1 \\ 0, & \boldsymbol{x}& \in &\Gamma\\ -\text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in & \Omega_2 \end{array} \right.

with \text{dist}(\boldsymbol{x}, \Gamma ) = \underset{\boldsymbol{y} \; \in \; \Gamma}{\min}( |\boldsymbol{x} - \boldsymbol{y}| ).

This function \phi had also the following property |\nabla\phi|=1.

Moreover, the unit normal vector \boldsymbol{n} outgoing from the interface and the curvature \mathcal{\kappa} can be obtained from the levelset function.

\boldsymbol{n}=\frac{\nabla\phi}{|\nabla\phi} \\ \mathcal{\kappa}=\nabla \cdot \boldsymbol{n}= \nabla \cdot \frac{\nabla\phi}{|\nabla\phi|}

Now we have exposed the levelset function, we need to define how the levelset will evolve and will spread into all the space. To do this, we use the following advection equation :

\partial_t\phi+\boldsymbol{u}\cdot\nabla\phi=0

where \boldsymbol{u} is an incompressible velocity field.

### Heaviside and Dirac functions

We define also the regularized Heaviside function H_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \dfrac{1}{2} \left(1+\dfrac{\phi}{\varepsilon}+\dfrac{\sin\left(\dfrac{\pi \phi}{\varepsilon}\right)}{\pi}\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 1, & \phi \geq \varepsilon \end{array} \right.

and the regularized Dirac function \delta_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \displaystyle\dfrac{1}{2 \varepsilon} \left(1+\cos\left(\dfrac{\pi \phi}{\varepsilon}\right)\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 0, & \phi \geq \varepsilon. \end{array} \right.

The first one gives a different value to each side of the interface ( here 0 in and 1 out ). The second one allow us to define quantities, with value different from 0 at the interface. A typical value of \varepsilon in literature is 1.5h where h is the mesh step size.

It should be noted that these functions allow us to determine respectively the volume and the surface of the interface by V^+_{\varepsilon} = \int_{\Omega} H_\varepsilon \\ S^{\Gamma}_{\varepsilon} = \int_{\Omega} \delta_\varepsilon

# Solid rotation of a slotted disk

We describe the benchmark proposed by Zalesak.

Computer codes, used for the acquisition of results, are from Vincent Doyeux.

## Problem Description

In order to test our interface propagation method, i.e. the levelset method \phi, we will study the rotation of a slotted disk into a square domain. The geometry can be represented as

Figure 1 : Initial Geometry.

We denote \Omega, the square domain [0,1]\times[0,1]. The center of the slotted disk is placed at (0.5,0.75).
To model the rotation, we will apply an angular velocity, centered in (0.5,0.5), as the disk is back to its initial position after t_f=628.

During this test, we observe three different errors to measure the quality of our method. With these values, two kinds of convergence will be studied : the time convergence, with different time step on an imposed grid and the space one, where the space discretization and the time step are linked by a relation. Several stabilization methods are used such as CIP ( Continuous Interior Penalty ) or SUPG ( Streamline-Upwing/Petrov-Galerkin ).

### Boundary conditions

We set a Neumann boundary condition on the boundary of the domain.

### Initial conditions

The velocity is imposed as \boldsymbol{u}=\left( \frac{\pi}{314} (50-y),\frac{\pi}{314} (x-50) \right)

Here is the velocity look in the square domain

Figure 2 : Imposed velocity in \Omega.

## Inputs

The following table displays the various fixed and variables parameters of this test-case.

 Name Description Nominal Value Units r disk radius 0.15 m l slot base 0.05 m h slot height 0.25 m t_f slotted disk rotation period 628 s

## Outputs

We observe during this benchmark three different errors.
First at all, the mass error, define by e_{\text{m}} = \frac{ \left| m_{\phi_f} - m_{\phi_0}\right| }{m_{\phi_0}} = \frac{ \left| \displaystyle \int_{\Omega} \chi( \phi_f < 0 ) - \displaystyle \int_{\Omega} \chi( \phi_0 < 0 ) \right| }{ \displaystyle \int_{\Omega} \chi( \phi_0 < 0 )}

where \chi is the characteristic function. \chi( f( \phi ) ) = \left\{ \begin{array}{rcl} 1 & \text{ if } & f( \phi ) \neq 0 \\ 0 & \text{ if } & f( \phi ) = 0 \end{array} \right.

However, mass is gain and loose at different emplacements on the mesh, and at the same time, with the level set method.

Secondly, the sign change error e_{\text{sc}} = \sqrt{ \int_\Omega \left( (1-H_0) - (1-H_f) \right)^2 }

with H_0=H_\epsilon(\phi_0) and H_f=H_\epsilon(\phi_f), H_\epsilon the smoothed Heaviside function of thickness 2ε.

This error is better to define the interface displacement. In fact, we can determine where \phi_0\phi_f<0, in other words where the interface has moved.

Finally, we define the classical L^2 error at the interface, as e_{L^2} = \sqrt{ \frac{1}{\displaystyle \int_\Omega \chi( \delta(\phi_0) > 0 ) } \int_\Omega (\phi_0 - \phi_f)^2 \chi( \delta(\phi_0) > 0 ) }.

## Discretization

### Time convergence

For this case, we set a fixed grid with mesh step size h=0.04, and so 72314 degree of freedom on a \mathbb{P}^1.

Then, after the disk made one round, we measure the errors obtained from two different discretizations ( BDF2 and Euler ) and compared them.

We repeat this with several time step dt\in \{2.14, 1, 0.5, 0.25, 0.20\}.

Only one stabilization method is used : SUPG

### Space convergence

We define the following relation, between time step and mesh step size : dt=C\frac{h}{U_{max}}

where C<1 constant and U_{max} the maximum velocity of \Omega.

From the definition of our velocity, U_{max} is reached at the farthest point from the center of \Omega. In this case, we have U_{max}=0.007, and we set C=0.8.

We use the BDF2 method for time discretization. As in time convergence, we wait one round of the disk to measure the errors and we repeat this test for these values of h: 0.32, 0.16, 0.08, 0.04.

We compare the results from different stabilization methods : CIP, SUPG, GLS ( Galerkin-Least-Squares ) and SGS ( Sub-Grid Scale ).

## Results

### Time convergence

 dt e_{L^2} e_{sc} e_m 2.14 0.0348851 0.202025 0.202025 1.00 0.0187567 0.147635 0.147635 0.5 0.0098661 0.10847 0.10847 0.25 0.008791 0.0782569 0.0782569 0.20 0.00803373 0.0670677 0.0670677
 dt e_{L^2} e_{sc} e_m 2.14 0.0118025 0.0906617 0.0492775 1.00 0.00436957 0.0445275 0.0163494 0.5 0.00173637 0.0216359 0.0100621 0.25 0.001003 0.0125971 0.00354644 0.20 0.000949343 0.0117449 0.00317368
Figure 2 : Mass error of Zalesac benchmark
Figure 3 : Sign change error of Zalesac benchmark
Figure 4 : L^2 error of Zalesac benchmark
Figure 5 : Slotted disk shape after a round

### Space convergence

 stab h e_{L^2} e_{sc} e_m CIP 0.32 0.0074 0.072 0.00029 0.16 0.0046 0.055 0.00202 0.08 0.0025 0.033 0.00049 0.04 0.0023 0.020 0.00110 SUPG 0.32 0.012 0.065 0.01632 0.16 0.008 0.049 0.07052 0.08 0.004 0.030 0.00073 0.04 0.001 0.018 0.00831 GLS 0.32 0.013 0.066 0.02499 0.16 0.008 0.051 0.05180 0.08 0.004 0.031 0.00805 0.04 0.001 0.019 0.00672 SGS 0.32 0.012 0.065 0.01103 0.16 0.008 0.050 0.07570 0.08 0.004 0.030 0.00084 0.04 0.001 0.018 0.00850
Figure 6 : Sign change error of Zalesac benchmark
Figure 7 : L^2 error of Zalesac benchmark
Figure 8 : Slotted disk shape after a round

### Conclusion

Let’s begin with time convergence results. Tables shows us that sign change error is better to define the quality of the chosen scheme than the mass error. In fact, the loss of mass somewhere can be nullify by a gain of mass elsewhere. Sign change error shows half an order gain from Euler scheme to BDF2 one, as L^2 errors show us a gain of one order. For the slotted disk shape, BDF2 uses the two previous iterations to obtain the current result, while Euler only need the previous iteration. This explain why we can see an asymmetrical tendency in the first one.

As for space convergence, the different stabilization methods we used give us the same convergence rate equals to 0.6, with close error values, for the sign change error. For the L^2 error case, it’s not as evident as the previous one. Aside the CIP stabilization method with a 0.6 convergence rate, the others show us a convergence rate of 0.9.

## Bibliography

References for this benchmark
• [Zalesak] Steven T. Zalesak, Fully multidimensional flux-corrected transport algorithms for fluids, Journal of Computational Physics, 1979.

• [Doyeux] Vincent Doyeux, Modelisation et simulation de systemes multi-fluides. Application aux ecoulements sanguins., Physique Numérique [physics.comp-ph], Université de Grenoble, 2014

:leveloffset:-1