Step by step into Feel++ Programming

Christophe Prud’homme <@prudhomm> v1.0, 2017/03/07

1. Output Directories

merge with 02-SettingUpEnvironment.adoc and rewrite pending!

Feel++ generates various files that are spread over various directories. For this tutorial, it would be beneficial to check the content of these files to familiarize yourself with Feel++.

1.1. Environment variables

Some of Feel++ behavior can be driven by environment variables such as

  • FEELPP_REPOSITORY

  • FEELPP_WORKDIR

Both variables should point to the same place. They define the root directory where the simulation results will be stored. By default they are set to $HOME/feel.

If you want to change the root directory where the results are stored, define e.g. FEELPP_REPOSITORY. For example in the docker image feel/apps:latest, it points to /feel.

For running,

docker run -it -v $HOME/feel:/feel feelpp/apps:latest

should get you this output

# Feel++ applications

This image provides the Feel++ applications in
 . advection
 . fluid mechanics
 . solid mechanics
 . fluid-structure interaction
 . thermodynamics
The testcases are in $HOME/Testcases/ and the results are in
$HOME/feel

Follow the steps described in

http://book.feelpp.org/content/Applications/Models/readme.html

to use the models.

feelpp@50381de2bd23:~$

and here is the result of echo $FEELPP_REPOSITORY in the docker image

feelpp@50381de2bd23:~pass:[\( echo \)]FEELPP_REPOSITORY
/feel

1.2. Results

  • Results $FEELPP_REPOSITORY/feel/<your_app_name>/np_1

1.3. Log files

Feel++ uses Google Glog.

  • Log files : $FEELPP_REPOSITORY/feel/<your_app_name>/np_1/logs

  • Mesh : $FEELPP_REPOSITORY/feel/<your_app_name>/np_1

  • Config files : $FEELPP_REPOSITORY/<your_build_folder>/doc/manual/tutorial

2. Setting up the Feel++ Environment

merge with 01-OutputDirectories.adoc and rewrite pending!

2.1. Minimal Example

Let’s begin with our first program using the Feel++ framework. To start, you include the Feel++ headers.

We use the C++ namespace to avoid Feel:: prefix before Feel++ objects.

We initialize the environment variables through the Feel++ Environment class, which can be found here.

#include <feel/feel.hpp>

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

    Environment env( _argc=argc, _argv=argv,
                     _about=about( _name="env",
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org") );
    std::cout << "proc " << Environment::rank()
              <<" of "<< Environment::numberOfProcessors()
              << std::endl;

}

and the config file

myapp-solver-type=cg
# myapp-pc-type=ilu

2.2. Adding options

We pass command line options using the Boost Program Options, library using the prefix po:: which is a Feel++ alias for the Boost::program_options namespace. To add a new Feel++ option, we must create a new Feel++ options_description. You must add the default Feel++ options and the new one that we choose here as a double value. Note that the default value will be assigned if not specified by the user.

2.3. Compilation execution and logs

To compile a tutorial, just use the GNU make command.

make feelpp_tut_<appname>

where <appname> is the name of the application you wish to compile (here, myapp). Go to the execution directory as specified in the program, and execute it.You can list the log files created :

ls /tmp/<your login>/feelpp/feelpp_tut_myapp/

If you open one of these log, you should be able to see your value and the processor number used to compute. You can run your application on several processors using MPI :

mpirun -np 2 feelpp_tut_myapp

Note that there will be one log for each processor in that case.

2.4. Config files

A config file can be parsed to the program to profile your options. The default config paths are,

  • current dir

  • $HOME/Feelpp/config/

  • $INSTALL_PREFIX/share/Feelpp/config/

then you have to write inside one of these folders a file called <app_name>.cfg or feelpp_<app_name>.cfg. For example, our myapp.cfg would look like :

value=0.53

Note that you can specify the config file through the option --config-file=<path>

It’s also possible to give several configuration files with the option --config-files <path1> <path2> <path3>

./feelpp_tut_myapp --config-files ex1.cfg ex2.cfg ex3.cfg

In the case where some options are duplicated in the files, the priority is given at the end :

  • ex3.cfg can overwrite options in ex2.cfg and ex1.cfg

  • ex2.cfg can overwrite options in ex1.cfg

All files in --config-files can overwrite options given by --config-file. And all options in the command line can overwrite all options given in cfg files.

2.5. Initializing PETSc, SLEPc and other third party libraries

PETSc is a suite of data structures and routines for the scalable (parallel) solution of scientific applications modeled by partial differential equations. It employs the MPI standard for parallelism.

Feel++ supports the PETSc framework, the Environment takes care of initializing the associated PETSc environment.

3. Loading a Mesh

The next step is to load a mesh.

The loadMesh function has a _name option set by default as the default value of the --gmsh.filename option that point either to a .geo, either to a .msh, or a .h5 file. Meshes in general are more detail into this section.

auto mesh=loadMesh( _mesh=new Mesh<Simplex<2>> );

3.1. Exporting the Mesh for visualisation

See this section for more details about exporting and visualizing meshes.

3.2. Implementation

An implementation reads as follows:

#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelvf/vf.hpp>

int main( int argc, char** argv )
{
    using namespace Feel;
    using namespace Feel::vf;
    // initialize Feel++ Environment
    Environment env( _argc=argc, _argv=argv,
                     _about=about( _name="mymesh" ,
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org" ) );
    // tag::mesh[]
    // create a mesh with GMSH using Feel++ geometry tool
    auto mesh = loadMesh(_mesh=new  Mesh<Simplex<2>>);
    // end::mesh[]

    LOG(INFO) << "mesh " << soption(_name="gmsh.filename") << " loaded";

    LOG(INFO) << "volume =" << integrate( _range=elements( mesh ), _expr=cst( 1. ) ).evaluate();
    LOG(INFO) << "surface = " << integrate( _range=boundaryfaces( mesh ), _expr=cst( 1. ) ).evaluate();
}

and the associated config file

[gmsh]
hsize=1e-1

4. Defining and using expressions

The next step is to construct a function space over the mesh.

4.1. Step by step explanations

  • We start by loading a Mesh in 2D.

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
  • then we define some expression through the command line of config file: g is a scalar field and f is a vector field, here is an example how to enter them :

./feelpp_tut_myexpression --a=3 --functions.g="a*x*y:x:y:a" --functions.f="{sin(pi*x),cos(pi*y)}:x:y"

You can print back the expression to the screen to check that everything is ok. You want to use as expression a*x+b*y, you have to define a and b as option (either in your code, either in the library).

  • then we compute the gradient of g and f.

    auto grad_g=grad<2>(g);
    auto grad_f=grad(f);
    std::cout << "grad(g)=" << grad_g << std::endl;
    std::cout << "grad(f)=" << grad_f << std::endl;
template argument are given to grad to specify the shape of the gradient: in the case of \(\nabla g\), it is \(1\times2\) and \(2\times 2\) for \(\nabla f\) since we are in 2D.
  • then we compute the laplacian of g and f.

    auto laplacian_g=laplacian(g);
    std::cout << "laplacian(g)=" << laplacian_g << std::endl;

    auto laplacian_f=laplacian(f);
    std::cout << "laplacian(f)=" << laplacian_f << std::endl;
  • then we compute the divergence of f.

    auto div_f=div(f);
    std::cout << "div(f)=" << div_f << std::endl;
  • and the curl of f

    auto curl_f=curl(f);
    std::cout << "curl(f)=" << curl_f << std::endl;
  • Finally we evaluate these expressions at one point given by the option x and y.

4.2. Implementation

An implementation reads as follows:

#include <feel/feelcore/environment.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feelvf/ginac.hpp>
using namespace Feel;

inline
po::options_description
makeOptions()
{
    po::options_description EXPRoptions( "DAR options" );
    EXPRoptions.add_options()
    ( "a", po::value<double>()->default_value( 1 ), "a parameter" )
    ( "b", po::value<double>()->default_value( 2 ), "a parameter" )
    ;
    return EXPRoptions;
}


int main(int argc, char**argv )
{
    Environment env( _argc=argc, _argv=argv,
                     _desc=makeOptions(),
                     _about=about(_name="myexpression",
                                  _author="Feel++ Consortium",
                                  _email="feelpp-devel@feelpp.org"));

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);

    auto g = expr(soption(_name="functions.g"));
    std::cout << "g=" << g << std::endl;

    auto f = expr<2,1>(soption(_name="functions.f"));
    std::cout << "f=" << f << std::endl;

    double aVal = doption("a")+doption("b");
    std::map<std::string,double> myMap; myMap["aVal"]=aVal;
    auto i = expr(soption("functions.i"),myMap);
    std::cout << "i=" << i << std::endl;

    auto grad_g=grad<2>(g);
    auto grad_f=grad(f);
    std::cout << "grad(g)=" << grad_g << std::endl;
    std::cout << "grad(f)=" << grad_f << std::endl;

    auto laplacian_g=laplacian(g);
    std::cout << "laplacian(g)=" << laplacian_g << std::endl;

    auto laplacian_f=laplacian(f);
    std::cout << "laplacian(f)=" << laplacian_f << std::endl;


    auto div_f=div(f);
    std::cout << "div(f)=" << div_f << std::endl;

    auto curl_f=curl(f);
    std::cout << "curl(f)=" << curl_f << std::endl;

    std::cout << "Evaluation  at  (" << doption("x") << "," << doption("y") << "):" << std::endl;
    std::cout << "           g(x,y)=" << g.evaluate() << std::endl;
    std::cout << "           f(x,y)=" << f.evaluate() << std::endl;
    std::cout << "           i(x,y)=" << i.evaluate() << std::endl;
    std::cout << "Gradient:\n";
    std::cout << "     grad(g)(x,y)=" << grad_g.evaluate() << std::endl;
    std::cout << "     grad(f)(x,y)=" << grad_f.evaluate() << std::endl;
    std::cout << "Divergence:\n";
    std::cout << "      div(f)(x,y)=" << div_f.evaluate() << std::endl;
    std::cout << "Curl:\n";
    std::cout << "     curl(f)(x,y)=" << curl_f.evaluate() << std::endl;
    std::cout << "Laplacian:\n";
    std::cout << "laplacian(g)(x,y)=" << laplacian_g.evaluate() << std::endl;
    std::cout << "laplacian(f)(x,y)=" << laplacian_f.evaluate() << std::endl;
}

and the associated config file

a=12
b=-1
[functions]
g=(a-x)*x+(a/b)*y^3:x:y:a:b
f={1,1}:x:y
i=(x-aVal)*y:x:y:aVal

4.3. Execution

$ ./feelpp_tut_myexpression

or

$ ./feelpp_tut_myexpression --a=3 --functions.g="<your_function>" --functions.f="<your_function>"

We start with the following function \(\)g=1\(\) and \(\)f=(1,1)\(\).

 $./feelpp_tut_myexpression --functions.g=1:x:y --functions.f="{1,1}:x:y"

and get something like this

g=1
f={1,1}
i=(x-aVal)*y
grad(g)=[[0,0]]
grad(f)=[[0,0],[0,0]]
laplacian(g)=[[0]]
laplacian(f)=[[0],[0]]
div(f)=[[0]]
curl(f)=[[0]]
Evaluation  at  (0,0):
           g(x,y)=1
           f(x,y)=1
1
           i(x,y)=-0
Gradient:
     grad(g)(x,y)=0 0
     grad(f)(x,y)=0 0
0 0
Divergence:
      div(f)(x,y)=0
Curl:
     curl(f)(x,y)=0
Laplacian:
laplacian(g)(x,y)=0
laplacian(f)(x,y)=0
0

The symbolic calculus system worked as expected.

5. Evaluating function

Once you have created an element, you may want to give it a value, that can depends on a lot of parameters ( mainly spaces, but others may apply ).

To do so, Feel++ relies on expressions. We may use various kind of expressions :

5.1. Built-in

Let’s begin with the evaluation of the expression \(sin(\pi x)\) on a unit circle.

First at all, we define the unit circle and its function space :

  // circle - geometrical order: 2
  auto mesh = unitCircle<2>();

  // circle - geometrical order: 1
  auto meshp1 = unitCircle<1>();

  auto Xh = Pch<2>( mesh );

Then the expression we would like to evaluate :

  auto myExpr = sin(pi*Px());

Px() refers to the variable \(\)x\(\) of our space.

With this,we can project it on our function space :

  auto v = project( _space=Xh, _range=elements(mesh),
      _expr=myExpr);

The expression will be evaluated on each point of our mesh.

In order to visualize the result, we create an exporter, named exhi, and add to it the projection.

  auto v = project( _space=Xh, _range=elements(mesh),
      _expr=myExpr);

5.1.1. Code

#include <feel/feel.hpp>

int main(int argc, char**argv )
{
  using namespace Feel;
  Environment env( _argc=argc, _argv=argv,
      _about=about(_name="myexporter",
        _author="Christophe Prud'homme",
        _email="christophe.prudhomme@feelpp.org"));

  // circle - geometrical order: 2
  auto mesh = unitCircle<2>();

  // circle - geometrical order: 1
  auto meshp1 = unitCircle<1>();

  auto Xh = Pch<2>( mesh );

  auto myExpr = sin(pi*Px());

  auto v = project( _space=Xh, _range=elements(mesh),
      _expr=myExpr);

  auto exhi = exporter( _mesh=mesh, _name="exhi" );
  auto exlo = exporter( _mesh=meshp1, _name="exlo" );
  auto exhilo = exporter( _mesh=lagrangeP1(_space=Xh)->mesh(),_name="exhilo");

  int max = 10; double dt = 0.1;
  double time = 0;
  for (int i = 0; i<max; i++)
  {
    exhilo->step( time )->add( "vhilo", v );
    exlo->step( time )->add( "vlo", v );
    exhi->step( time )->add( "vhi", v );
    time += dt;
    exhi->save();
    exlo->save();
    exhilo->save();
  }
}

The list of the Feel++ Keyword is here.

5.2. Hard Coded

In this second method, we will use Functor :

struct Functor
{
    static const size_type context = vm::JACOBIAN|vm::POINT;
    typedef double value_type;
    typedef Feel::uint16_type uint16_type;
    static const uint16_type rank = 0;
    static const uint16_type imorder = 1;
    static const bool imIsPoly = true;
    myFunctor( int i ): coord(i){}
    double operator()( uint16_type a, uint16_type b, ublas::vector<double> const& x, ublas::vector<double> const& n ) const
        {
            return xtag::coord[];
        }
    int coord = 0;
};

We create a unit square meshed by triangles and we define the associated function space :

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    auto Xh= Pch<2>(mesh);

From this space, we can define two elements, here one equals to the variable \(x\) and the other to the variable \(y\), obtain from Functor class.

    auto u = Xh.element(idf(Functor(0))); // Will contain x
    auto v = Xh.element(idf(Functor(1))); // Will contain y

The data exportation is the final step to visualize our expression \(x\) and \(y\) on the defined mesh.

    auto ex1 = exporter(_mesh=mesh,_name="ex1");
    ex1->add("u",u);
    ex1->add("v",v);
    ex1->save();

5.2.1. Code

The complete code reads as follows

#include <feel/feelcore/environment.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feeldiscr/pch.hpp>
#include <feel/feeldiscr/pchv.hpp>
#include <feel/feelvf/projectors.hpp>
#include <feel/feeldiscr/operatorlagrangep1.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelvf/function.hpp>
using namespace Feel;

namespace Feel {
struct Functor
{
    static const size_type context = vm::JACOBIAN|vm::POINT;
    typedef double value_type;
    typedef Feel::uint16_type uint16_type;
    static const uint16_type rank = 0;
    static const uint16_type imorder = 1;
    static const bool imIsPoly = true;
    myFunctor( int i ): coord(i){}
    double operator()( uint16_type a, uint16_type b, ublas::vector<double> const& x, ublas::vector<double> const& n ) const
        {
            return xtag::coord[];
        }
    int coord = 0;
};
}

int main(int argc, char**argv )
{
    Environment env( _argc=argc, _argv=argv,
                     _about=about(_name="myfunctor",
                                  _author="Feel++ Consortium",
                                  _email="feelpp-devel@feelpp.org"));


    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    auto Xh= Pch<2>(mesh);


    auto u = Xh.element(idf(Functor(0))); // Will contain x
    auto v = Xh.element(idf(Functor(1))); // Will contain y

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

The list of the Feel++ Keyword is here.

6. Visualizing functions over a mesh

The next step is to visualize function over the mesh. The source code that generate our output is available in myexporter.cpp, presented in the previous section.

You can visualize data via :

  • link::https://www.ceisoftware.com/[ensight]

  • link::http://www.paraview.org/[paraview]

  • link::http://geuz.org/gmsh[gmsh]

The results files are in

  • $HOME/feel/myexporter/np_1

or

  • $FEELPP_WORKDIR/feel/myexporter/np_1.

We discriminate output directories based on the name of the simulation (the parameter _name in the environment), the number of process (mpirun -np …​) and the type of the chosen exporter

--exporter.format={ensight|ensightgold|gmsh|...}

7. Spaces and elements

You’ve learned how to discretize the space you want to compute on. You now have to learn how to define and use function spaces and elements of functions spaces. For advanced informations on this subject, you can look in the Function Space documentation.

7.1. Constructing a function space

  • Loading a Mesh in 2D

auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
  • For basic function spaces, we have predetermined constructors:

auto Xh = Pch<2>( mesh );
  • Defining an element

auto u = Xh->element( "u" );
auto w = Xh->element( "w" );

One can also use :

  • Pdh<ORDER>(mesh) : Polynomial Discontinuous

  • Pvh<ORDER>(mesh) : Polynomial Continuous Vectorial

  • Pdhv<ORDER>(mesh) : Polynomial Discontinuous Vectorial

  • Pchm<ORDER>(mesh) : Polynomial Continuous Matrix

  • Ned1h<ORDER>(mesh) : Nedelec function spaces

7.2. Implementation

The implementation reads are follows

#include <feel/feel.hpp>
using namespace Feel;

int main( int argc, char** argv )
{

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

    // create the mesh
    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);

    // function space \( X_h \) using order 2 Lagrange basis functions
    auto Xh = Pch<2>( mesh );

    auto g = expr<4>( soption(_name="functions.g"));
    auto gradg = grad<3>(g);

    // elements of \( u,w \in X_h \)
    auto u = Xh->element( "u" );
    auto w = Xh->element( "w" );

		// build the interpolant of u
    u.on( _range=elements( mesh ), _expr=g );

		// build the interpolant of the interpolation error
    w.on( _range=elements( mesh ), _expr=idv( u )-g );

    // compute L2 norms \(||\cdot||_{L^2}\)
    double L2g = normL2( elements( mesh ), g );
    double H1g = normL2( elements( mesh ), _expr=g,_grad_expr=gradg );
    double L2uerror = normL2( elements( mesh ), ( idv( u )-g ) );
    double H1uerror = normH1( elements( mesh ), _expr=( idv( u )-g ),
                              _grad_expr=( gradv( u )-gradg ) );
    if ( Environment::isMasterRank() )
    {
        std::cout << "||u-g||_0 = " << L2uerror/L2g << std::endl;
        std::cout << "||u-g||_1 = " << H1uerror/H1g << std::endl;
    }

    // export for post-processing
    auto e = exporter( _mesh=mesh );

	 	// save interpolant
    e->add( "g", u );

	 	// save interpolant of interpolation error
    e->add( "u-g", w );

    e->save();

}

and the associated config file

[gmsh]
# default is hypercube
shape=hypercube
# you can try ellipsoid
#shape=ellipsoid

[functions]
g=x*y*(y*2.2+4.3)*(2.1*x-1.3):x:y

8. Computing integrals over mesh

The next step is to compute integrals over the mesh ( See this for detailed methods ).

8.1. Step by step explanations

  • We start by loading a Mesh in 2D

  • then we define the Feel++ expression that we are going to integrate using the soption function that retrieves the command line option string functions.g. We then transform this string into a Feel++ expression using expr().

  • then we compute two integrals over the domain and its boundary respectively

    • \(\int_\Omega g\)

    • \(\int_{\partial \Omega} g\)

  • and we print the results to the screen.

Only the rank 0 process (thanks to Environment) isMasterRank() prints to the screen as the result is the same over all mpi processes if the application was run in parallel. Note also that the code actually prints the expression passed by the user through the command line option functions.g.

8.2. Some results

We start with the following function \(g=1\). Recall that by default the domain is the unit square, hence the \(\int_\Omega g\) and \(\int_{\partial\Omega} g\) should be equal to 1 and 4 respectively.

./feelpp_tut_myintegrals --functions.g=1
int_Omega 1 = 1
int_{boundary of Omega} 1 = 4

Now we try \(g=x\). We need to tell Feel++ what are the symbols associated with the expression: here the symbol is x and it works as follows

./feelpp_tut_myintegrals --functions.g=x:x
int_Omega x = 0.5
int_{boundary of Omega} x = 2
remember that there is a separator : between the expression and each symbol composing it.

Now we try \(g=x y\). We need to tell Feel++ what are the symbols associated with the expression: here the symbol is x and y and it works as follows

./feelpp_tut_myintegrals --functions.g=x*y:x:y
int_Omega x*y = 0.25
int_{boundary of Omega} x*y = 1

More complicated functions are of course doable, such as \(g=\sin( x y ).\)

./feelpp_tut_myintegrals --functions.g="sin(x*y):x:y"
int_Omega sin(x*y) = 0.239812
int_{boundary of Omega} sin(x*y) = 0.919395

Here is the last example in parallel over 4 processors which returns, of course, the exact same results as in sequential

mpirun -np 4 ./feelpp_doc_myintegrals --functions.g="sin(x*y):x:y"
int_Omega sin(x*y) = 0.239812
int_{boundary of Omega} sin(x*y) = 0.919395

Finally we can change the type of domain and compute the area and perimeter of the unit disk as follows

./feelpp_doc_myintegrals --functions.g="1:x:y" --gmsh.domain.shape=ellipsoid --gmsh.hsize=0.05
int_Omega 1 = 0.784137
int_{boundary of Omega} 1 = 3.14033

Note that we don’t get the exact results due to the fact that [stem]:[\Omega_h = \cup_{K \in \mathcal{T}_h} K] which we use for the numerical integration is different from the exact domain \(\Omega = \{ (x,y)\in \mathbb{R}^2 | x^2+y^2 < 1\}\).

8.3. Implementation

To compile just type

$ ./feelpp_tut_myintegrals

The complete code reads as follows :

#include <feel/feel.hpp>


using namespace Feel;

int
main( int argc, char** argv )
{
    // Initialize Feel++ Environment
    Environment env( _argc=argc, _argv=argv,
                     _about=about( _name="myintegrals" ,
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org" ) );

    /// [mesh]
    // create the mesh (specify the dimension of geometric entity)
    auto mesh = loadMesh( _mesh=new Mesh<Simplex<2>> );
    /// [mesh]

    /// [expression]
    // our function to integrate
    auto g = expr( soption(_name="functions.g") );
    /// [expression]

    /// [integrals]
    // compute integral of g (global contribution): \(\int_{\partial \Omega} g\)
    auto intf_1 = integrate( _range = elements( mesh ),
                                 _expr = g ).evaluate();

    // compute integral g on boundary: \( \int_{\partial \Omega} g \)
    auto intf_2 = integrate( _range = boundaryfaces( mesh ),
                             _expr = g ).evaluate();

    // compute integral of grad f (global contribution): \( \int_{\Omega} \nabla g \)
    auto grad_g = grad<2>(g);
    auto intgrad_f = integrate( _range = elements( mesh ),
                                _expr = grad_g ).evaluate();

    // only the process with rank 0 prints to the screen to avoid clutter
    if ( Environment::isMasterRank() )
        std::cout << "int_Omega " << g << " = " << intf_1  << std::endl
                  << "int_{boundary of Omega} " << g << " = " << intf_2 << std::endl
                  << "int_Omega grad " << g << " = "
                  << "int_Omega  " << grad_g << " = "
                  << intgrad_f  << std::endl;
    /// [integrals]
}
/// [all]
[functions]
g=x*y*(y*2.2+4.3)*(2.1*x-1.3):x:y

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

9.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\]

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

10. Using a backend

10.1. Introduction

After the discretization process, one may have to solve a (non) linear system. Feel++ interfaces with PETSc/SLEPc and Eigen3. Consider this system

\[A x = b\]

We call Backend an object that manages the solution strategy to solve it. Some explanation are available at Solver and Preconditioner.

Feel++ provides a default backend that is mostly hidden to the final user. In many examples, you do not have to take care of the backend. You change the backend behavior via the command line or config files. For example

./feelpp_doc_mybackend --backend.pc-type=id

will use the identity matrix as a right preconditionner for the default backend. The size of the preconditionner will be defined from the size of the A matrix.

If you try to solve a different system \(A_1 y= c\) (in size) with the same backend or the default without rebuilding it, it will fail.

backend(_rebuild=true)->solve(_matrix=A1,_rhs=c,_sol=y);

Each of that options can be retrieved via the --help-lib argument in the command line.

10.1.1. Non default Backend

You may need to manage more than one backend in an application: you have different systems to solve and you want to keep some already computed objects such as preconditioners.

  • The default backend is in fact an unnamed backend: in order to distinguish between backend you have to name them. For example

	po::options_description app_options( "MyBackend options" );
	app_options.add(backend_options("myBackend"));

	Environment env(_argc=argc, _argv=argv,
	                _desc = app_options,
	                _about=about(_name="mybackend",
	                             _author="Feel++ Consortium",
	                             _email="feelpp-devel@feelpp.org"));
  • After that, you create the backend object:

	// create a backend
	boost::shared_ptr<Backend<double>> myBackend(backend(_name="myBackend"));
the backend’s name has to match the name you gave at the options step.
  • Then, you load meshes, creates spaces etc. At solve time, or you solve with the default backend:

	// solve a(u,v)=l(v)
	if ( Environment::isMasterRank() )
		std::cout << "With default backend\n";
	a.solve(_rhs=l,_solution=u1); // Compute with default backend

One of the important backend option is to be able to monitor the residuals and iteration count

./feelpp_tut_mybackend --pc-type=id --ksp-monitor=true --myBackend.ksp-monitor=true
  • Finally you can create a named backend:

	// solve a(u,v)=l(v)
	if ( Environment::isMasterRank() )
		std::cout << "With named backend\n";
	a.solveb(_rhs=l,_solution=u2, _backend=myBackend); // Compute with myBackend

10.2. Implementation

The implementation reads as follows:

#include <feel/feel.hpp>
using namespace Feel;

int main(int argc, char**argv )
{
	po::options_description app_options( "MyBackend options" );
	app_options.add(backend_options("myBackend"));

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

	// create a backend
	boost::shared_ptr<Backend<double>> myBackend(backend(_name="myBackend"));

	// create the mesh
	auto mesh = loadMesh(_mesh=new Mesh<Simplex< 2 > > );

	// function space
	auto Vh = Pch<2>( mesh );

	// element in Vh
	auto u  = Vh->element();
	auto u1 = Vh->element();
	auto u2 = Vh->element();

	// left hand side
	auto a = form2( _trial=Vh, _test=Vh );
	a = integrate(_range=elements(mesh),
	              _expr=expr(soption("functions.alpha"))*trace(gradt(u)*trans(grad(u))) );

	// right hand side
	auto l = form1( _test=Vh );
	l = integrate(_range=elements(mesh),
	              _expr=expr(soption("functions.f"))*id(u));

	// BC
	a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
	      _expr=expr(soption("functions.g")));

	// solve a(u,v)=l(v)
	if ( Environment::isMasterRank() )
		std::cout << "With default backend\n";
	a.solve(_rhs=l,_solution=u1); // Compute with default backend
	// solve a(u,v)=l(v)
	if ( Environment::isMasterRank() )
		std::cout << "With named backend\n";
	a.solveb(_rhs=l,_solution=u2, _backend=myBackend); // Compute with myBackend

	// save results
	auto e = exporter( _mesh=mesh );
	e->step(0) -> add( "u", u1 );
	e->step(1) -> add( "u", u2 );
	e->save();
}

and the associated config file

ksp-monitor=true
ksp-rtol=1e-5
backend.verbose=true
pc-type=id

[functions]
alpha=1
f=x+y:x:y
g=sin(x):x

[myBackend]
backend.verbose=true
ksp-monitor=true
ksp-rtol=1e-5
pc-type=lu

11. Defining a Model

11.1. Introduction

It is well known an equation can describe a huge range of physical problems. Each of theses problems will have a particular environment, but the equation to solve will be the same. To make our program applicable to theses range of problem, we have defined a model. Models definitions can be retrieve in this section.

11.2. What is a model

A model is defined by :

  • a Name

  • a Description

  • a Model

  • Parameters

  • Materials

  • Boundary Conditions

  • Post Processing

11.2.1. Parameters

A parameter is a non physical property for a model.

11.2.2. Materials

To retrieve the materials properties, we use :

ModelMaterials materials = model.materials();

11.2.3. BoundaryConditions

Thanks to GiNaC, we handle boundary conditions (Dirichlet, Neumann, Robin) as expression. You have to indicate in the json file the quantity to handle (velocity, pressure…​) and the associated expression.

map_scalar_field<2> bc_u { model.boundaryConditions().getScalarFields<2>("heat","dirichlet") };

We can apply theses boundary condition this way

  for(auto it : bc_u){
    if(boption("myVerbose") && Environment::isMasterRank() )
      std::cout << "[BC] - Applying " << it.second << " on " << it.first << std::endl;
    a+=on(_range=markedfaces(mesh,it.first), _rhs=l, _element=u, _expr=it.second );
  }

11.2.4. Code

  for(auto it : materials)
  {
    auto mat = material(it);
    if(boption("myVerbose") && Environment::isMasterRank() )
      std::cout << "[Materials] - Laoding data for " << it.second.name() << " that apply on marker " << it.first  << " with diffusion coef ["
#if MODEL_DIM == 3
        << "[" << it.second.k11() << "," << it.second.k12() << "," << it.second.k13() << "],"
        << "[" << it.second.k12() << "," << it.second.k22() << "," << it.second.k23() << "],"
        << "[" << it.second.k13() << "," << it.second.k23() << "," << it.second.k33() << "]]"
#else
        << "[" << it.second.k11() << "," << it.second.k12() << "],"
        << "[" << it.second.k12() << "," << it.second.k22() << "]]"
#endif
        << std::endl;
    k11.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k11()));
    k12.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k12()));
    k22.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k22()));
#if MODEL_DIM == 3
    k13 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k13());
    k23 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k23());
    k33 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k33());
#endif
  }
#if MODEL_DIM == 2
  a += integrate(_range=elements(mesh),_expr=inner(mat<MODEL_DIM,MODEL_DIM>(idv(k11), idv(k12), idv(k12), idv(k22) )*trans(gradt(u)),trans(grad(v))) );
#else
  a += integrate(_range=elements(mesh),_expr=inner(mat<MODEL_DIM,MODEL_DIM>(idv(k11), idv(k12), idv(k13), idv(k12), idv(k22), idv(k23), idv(k31), idv(k32), idv(k33))*trans(gradt(u)),trans(grad(v))) );
#endif

11.2.5. PostProcessing

explanation pending.

11.2.6. Example

We have set up an example : an anisotropic laplacian.

#include <feel/feel.hpp>
#include <feel/feelmodels/modelproperties.hpp>
int main(int argc, char**argv )
{
  using namespace Feel;
  po::options_description laplacianoptions( "Laplacian options" );
  laplacianoptions.add_options()
    ("myVerbose", po::value< bool >()-> default_value( true ), "Display information during execution")
    ;
  Environment env( _argc=argc, _argv=argv,
      _desc=laplacianoptions,
      _about=about(_name="aniso_laplacian",
        _author="Feel++ Consortium",
        _email="feelpp-devel@feelpp.org"));
  ModelProperties model; // Will load --mod-file
  map_scalar_field<2> bc_u { model.boundaryConditions().getScalarFields<2>("heat","dirichlet") };
  ModelMaterials materials = model.materials();
  if(boption("myVerbose") && Environment::isMasterRank() )
    std::cout << "Model " << Environment::expand( soption("mod-file")) << " loaded." << std::endl;
  auto f = expr( soption(_name="functions.f"), "f" );
  auto mesh = loadMesh(_mesh=new Mesh<Simplex<MODEL_DIM>>);
  auto Vh = Pch<2>( mesh );
  auto u = Vh->element();
  auto v = Vh->element();
  auto k11 = Vh->element();
  auto k12 = Vh->element();
  auto k22 = Vh->element();
#if MODEL_DIM == 3
  auto k13 = Vh->element();
  auto k23 = Vh->element();
  auto k33 = Vh->element();
#endif
  auto a = form2( _trial=Vh, _test=Vh);
  auto l = form1( _test=Vh );
  l = integrate(_range=elements(mesh),_expr=f*id(v));
  for(auto it : materials)
  {
    auto mat = material(it);
    if(boption("myVerbose") && Environment::isMasterRank() )
      std::cout << "[Materials] - Laoding data for " << it.second.name() << " that apply on marker " << it.first  << " with diffusion coef ["
#if MODEL_DIM == 3
        << "[" << it.second.k11() << "," << it.second.k12() << "," << it.second.k13() << "],"
        << "[" << it.second.k12() << "," << it.second.k22() << "," << it.second.k23() << "],"
        << "[" << it.second.k13() << "," << it.second.k23() << "," << it.second.k33() << "]]"
#else
        << "[" << it.second.k11() << "," << it.second.k12() << "],"
        << "[" << it.second.k12() << "," << it.second.k22() << "]]"
#endif
        << std::endl;
    k11.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k11()));
    k12.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k12()));
    k22.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k22()));
#if MODEL_DIM == 3
    k13 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k13());
    k23 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k23());
    k33 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k33());
#endif
  }
#if MODEL_DIM == 2
  a += integrate(_range=elements(mesh),_expr=inner(mat<MODEL_DIM,MODEL_DIM>(idv(k11), idv(k12), idv(k12), idv(k22) )*trans(gradt(u)),trans(grad(v))) );
#else
  a += integrate(_range=elements(mesh),_expr=inner(mat<MODEL_DIM,MODEL_DIM>(idv(k11), idv(k12), idv(k13), idv(k12), idv(k22), idv(k23), idv(k31), idv(k32), idv(k33))*trans(gradt(u)),trans(grad(v))) );
#endif
  for(auto it : bc_u){
    if(boption("myVerbose") && Environment::isMasterRank() )
      std::cout << "[BC] - Applying " << it.second << " on " << it.first << std::endl;
    a+=on(_range=markedfaces(mesh,it.first), _rhs=l, _element=u, _expr=it.second );
  }
  a.solve(_rhs=l,_solution=u);
  auto e = exporter( _mesh=mesh );
  for(int i = 0; i < 3; i ++){
    for(auto const &it : model.postProcess()["Fields"] )
    {
      if(it == "diffused")
        e->step(i)->add("diffused",u);
      else if(it == "k11")
        e->step(i)->add("k11",k11);
      else if(it == "k12")
        e->step(i)->add("k12",k12);
      else if(it == "k11")
        e->step(i)->add("k22",k22);
#if MODEL_DIM == 3
      else if(it == "k13")
        e->step(i)->add("k13",k13);
      else if(it == "k11")
        e->step(i)->add("k23",k23);
      else if(it == "k33")
        e->step(i)->add("k33",k33);
#endif
    }
    e->save();
  }
  return 0;
}

Programming in Feel++

Unresolved directive in /var/lib/buildkite-agent/builds/sd-87660-1/feelpp/book-dot-feelpp-dot-org/docs/man/06-programming/README.adoc - include::https://raw.githubusercontent.com/feelpp/feelpp/develop/CONTRIBUTING.adoc[]

1. Feel++ Coding Styles

Christophe Prud’homme <@prudhomm> v1.0, 2017/03/07

This is an overview of the coding conventions we use when writing Feel++ code.

1.1. Clang Format

clang-format is a powerful tool to reformat your code according to rules defined in a .clang-format file at the toplevel directory of your software.

Feel++ has such a file and define the indentation, space and breaks rules defined later on.

For clang-format to function properly, follow the Comments rules.

To apply Feel++ rules on a file a.cpp in a Feel++ sub-directory, type

clang-format a.cpp

to dump the results of the reformating to the standard output or type

clang-format -i a.cpp

which will replace a.cpp by the reformated file.

be careful when reformating, make sure nobody is working on that file, to avoid creating possibly massive conflicts with the persons currently modifying the same code when they get merged.

1.2. Header files and multiple inclusions

To avoid multiple inclusions, wrap every header files using the following technique

// say we have myheader.hpp
#if !defined(FEELPP_MYHEADER_HPP)
#define FEELPP_MYHEADER_HPP 1

// your header here...

#endif // FEELPP_MYHEADER_HPP

more details here

1.3. Naming Convention

In Feel++, we basically follow the same naming conventions as in Qt and KDE.

Class names starts with a capital. The rest is in camel case. Function names starts with a lower case, but the first letter of each successive word is capitalized.

Functions and classes should be in the Feel namespace.

The prefix set is used for setters, but the prefix get is not used for accessors. Accessors are simply named with the name of the property they access. The exception is for accessors of a boolean which may start with the prefix is.

Acronyms are lowercased too. Example: Url instead of URL and isFemEnabled() instead of isFEMEnabled()

Accessors should usually be const.

This example shows some possible functions names

class A
{
public:
    void setMatrix(const Matrix& c);
    Matrix matrix() const;
    void setDiagonal(bool b);
    bool isDiagonal() const;
};

1.4. Reserved identifiers

1.5. Indentation

  • 4 spaces are used for indentation but not in namespace

  • Spaces, not tabs!

  • Suggestion: use emacs and [http://emacswiki.org/emacs/dirvars.el dirvars.el], here is the content of .emacs-dirvars in top Feel++ directory

    indent-tabs-mode: nil
    tab-width: 4
    c-basic-offset: 4
    evaluate: (c-set-offset 'innamespace '0)
    show-trailing-whitespace: t
    indicate-empty-lines: t
    evaluate: (add-hook 'write-file-hooks 'delete-trailing-whitespace)
namespace Feel
{
// no space indentation in namespace
Class A
{
    // 4 spaces indentation
    A() {};
    void f();
};
}

1.6. Comments

Use C++ style comment // rather than /* */. It uses less characters but also it is easier to reflow using clang format.


/* Wrong
  the doc
  */

// Correct
// the doc

1.6.1. Doxygen

Doxygen is the tool to document Feel++ and create a reference manual

Use //! to comment function, variables, classes rather than /** */, it allows to reflow comments using clang format.

//!
//! @brief the class
//! @author me <me@email>
//!
class TheClass
{
public:
  //! constructor
  TheClass() {}

private:
  //! member
  int member;
};

//! the function
void thefunction() {}
Feel++ used to promote /** */ but this is no longer the case. The comment style will be updated progressively to match the new style using //!

1.7. Declaring variables

  • Declare each variable on a separate line

  • Avoid short (e.g. a, rbarr, nughdeget) names whenever possible

  • Single character variable names are only okay for counters and temporaries, where the purpose of the variable is obvious

  • Wait when declaring a variable until it is needed

// Wrong
int a, b;
char __c, __d;

// Correct
int height;
int width;
char __nameOfThis;
char __nameOfThat;
  • Variables and functions start with a lower-case letter. Each consecutive word in a variable’s or function’s name starts with an upper-case letter

  • Avoid abbreviations


// Wrong
short Cntr; char ITEM_DELIM = '';

// Correct
short counter; char itemDelimiter = '';
  • Classes always start with an upper-case letter.


// Wrong
class meshAdaptation {};

// Correct
class MeshAdaptation {};
  • Non-static data members name of structures and classes always start with M_ . M stands for Member. The rational behind this is for example :

    • to be able to immediately see that the data is a member of a class or a struct

    • to easily search and query-replace


// Wrong
class meshAdaptation { std::vector directions_; };

// Correct
class MeshAdaptation { std::vector M_directions; };
  • Static data members name of structures and classes always start with S_ . S stands for Static. The rational behind this is for example :

    • to be able to immediately see that the data is a static member of a class or a struct

    • to easily search and query-replace


// Wrong
class meshAdaptation { static std::vector directions_; };

// Correct
class MeshAdaptation { static std::vector S_directions; };

1.8. Whitespace

  • Use blank lines to group statements together where suited

  • Always use only one blank line

  • Always use a single space after a keyword and before a curly brace.


// Correct
if (foo) { }

// Wrong
if(foo) { }
  • For pointers or references, always use a single space between the type and or &, but no space between the or & and the variable name.


char *x;
const std::string &myString;
const char * const y = "hello";
  • Surround binary operators with spaces.

  • No space after a cast.

  • Avoid C-style casts when possible.


// Wrong
char* blockOfMemory = (char* ) malloc(data.size());

// Correct
char *blockOfMemory = reinterpret_cast(malloc(data.size()));

1.9. Braces

  • As a base rule, the left curly brace goes on the same line as the start of the statement:


// Wrong
if (codec) { }

// Correct
if (codec) { }
  • Function implementations and class declarations always have the left brace on the start of a line:


static void foo(int g) { std::cout << g << "" }

class Moo { };
  • Use curly braces when the body of a conditional statement contains more than one line, and also if a single line statement is somewhat complex.

// Wrong
if (address.isEmpty()) { return false; }

for (int i = 0; i < 10; ++i) { std::cout << "i=" << i << ""; }

// Correct
if (address.isEmpty()) return false;

for (int i = 0; i < 10; ++i) std::cout << "=" << i << "";
  • Exception 1: Use braces also if the parent statement covers several lines / wraps

// Correct
if (address.isEmpty() || !isValid() || !codec)
{
    return false;
}
  • Exception 2: Use braces also in if-then-else blocks where either the if-code or the else-code covers several lines

// Wrong
if (address.isEmpty())
  return false;
else
{
  std::cout << address << ""; ++it;
 }

// Correct
if (address.isEmpty())
{
    return false;
}
else
{
    std::cout << address << ""; ++it;
}

// Wrong
if (a) if (b) ... else ...

// Correct
if (a) { if (b) ... else ... }
  • Use curly braces when the body of a conditional statement is empty


// Wrong
while (a);

// Correct
while (a) {}

1.10. Parentheses

  • Use parentheses to group expressions:


// Wrong
if (a && b || c)

// Correct
if ((a && b) || c)

// Wrong
a + b & c

// Correct
(a + b) & c

1.11. Switch statements

  • The case labels are in the same column as the switch

  • Every case must have a break (or return) statement at the end or a comment to indicate that there’s intentionally no break, unless another case follows immediately.

switch (myEnum)
{
  case Value1:
    doSomething();
    break;
  case Value2:
  case Value3:
    doSomethingElse(); // fall through
  default:
    defaultHandling();
    break;
}

1.12. Line breaks

  • Keep lines shorter than 100 characters; insert breaks if necessary.

  • Commas go at the end of a broken line; operators start at the beginning of the new line. An operator at the end of the line is easy to not see if your editor is too narrow.

// Correct
if (longExpression + otherLongExpression + otherOtherLongExpression) { }

// Wrong
if (longExpression + otherLongExpression + otherOtherLongExpression) { }

1.13. Inheritance and the virtual keyword

When reimplementing a virtual method, do not put the virtual keyword in the header file.