Solve a partial differential equation

With all the previously notions we approach, the definition of a partial differential equation and boundary conditions are our next step. More details on these aspects can be retrieve at this page.

1. Variational formulation

This example refers to a laplacian problem, define by

Strong formulation
\[-\Delta u=1 \text{ in } \Omega=[0,1]^2, \quad u=1 \text{ on } \partial \Omega\]

After turning the Strong formulation into its weak form, we have

\[\int_\Omega \nabla u \cdot \nabla v = \int_\Omega v,\quad u=1 \text{ on } \partial \Omega\]

where \(u\) is the unknown and \(v\) a test function. The left side is known as the bilinear form \(a\) and the right side is the linear form \(l\).

\[a(u,v) = l(v), \quad u=1 \text{ on } \partial \Omega\]

2. Implementation

The steps to implement this problem are

  • Loading a 2D mesh, creating the function space \(V_h\), composed of piecewise polynomial functions of order 2, and its associated elements

 auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
 auto Vh = Pch<2>( mesh );
 auto u = Vh->element();
 auto v = Vh->element();

  • Define the linear form \(l\) with test function space \(V_h\)

auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
              _expr=id(v));

  • Define the bilinear form \(a\) with \(V_h\) as test and trial function spaces

auto a = form2( _trial=Vh, _test=Vh);
a = integrate(_range=elements(mesh),
              _expr=gradt(u)*trans(grad(v)) );

form1 and form2 are used to define respectively the left and right side of our partial differential equation.

  • Add Dirichlet boundary condition on u

a+=on(_range=boundaryfaces(mesh),
      _rhs=l, _element=u, _expr=cst(0.) );

We impose, in this case, \(u=0\) on \(\partial\Omega\), with the keyword on.

  • Solving the problem

a.solve(_rhs=l,_solution=u);

  • Exporting the solution

auto e = exporter( _mesh=mesh );
e->add( "u", u );
e->save();

The complete code reads as follows :

#include <feel/feel.hpp>

int main(int argc, char**argv )
{
     using namespace Feel;

     Environment env( _argc=argc, _argv=argv,
                      _about=about(_name="qs_laplacian",
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org"));

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    auto Vh = Pch<2>( mesh );
    auto u = Vh->element();
    auto v = Vh->element();

    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=id(v));

    auto a = form2( _trial=Vh, _test=Vh);
    a = integrate(_range=elements(mesh),
                  _expr=gradt(u)*trans(grad(v)) );
    a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=cst(0.) );
    a.solve(_rhs=l,_solution=u);

    auto e = exporter( _mesh=mesh );
    e->add( "u", u );
    e->save();
}

and the corresponding config file

[gmsh]
hsize=1e-1