[programming]

Solving for the Laplacian

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

Table of Contents

The Laplacian

Feel++ Book Contributors https://github.com/feelpp/book.feelpp.org/graphs/contributors

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)\).

The weak formulation reads:

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

    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),
                  _expr=mu*inner(gradt(u),grad(v)) );
    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();
            double h1 = normH1(_range=elements(mesh), _expr=idv(u)-expr(solution), _grad_expr=gradv(u)-grad<2>(expr(solution)) );
            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 );
    e->addRegions();
    e->add( "uh", u );
    if ( checker().check() )
    {
        v.on(_range=elements(mesh), _expr=solution );
        e->add( "solution", v );
    }
    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.

ucircle

meshCircle

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

ufeelpp2d

/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

ufeelpp3d

meshfeelpp3d

Solution \(u_\delta\)

Mesh

6. All Generic Programming Examples

[programming] · Solving for the Laplacian