1. Computing integrals over mesh

The next step is to compute integrals over the mesh.

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

1.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 \(\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\}\).

1.3. Implementation

To compile just type

$ ./feelpp_tut_myintegrals

The complete code reads as follows :

#include <feel/feelcore/environment.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feelvf/integrate.hpp>
#include <feel/feelvf/ginac.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
    Feel::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