## 1. Computing integrals over mesh

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

### 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 [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\}$.

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