# Integration

You should be able to create a mesh now. If it is not the case, get back to the section Mesh.

Prerequisites

## 1. Integrals

Feel++ provide the integrate() function to define integral expressions which can be used to compute integrals, define linear and bi-linear forms.

### 1.1. Interface

``  integrate( _range, _expr, _quad, _geomap );``

Please notice that the order of the parameter is not important, these are `boost` parameters, so you can enter them in the order you want. To make it clear, there are two required parameters and 2 optional and they of course can be entered in any order provided you give the parameter name. If you don’t provide the parameter name (that is to say `_range` = or the others) they must be entered in the order they are described below.

Required parameters:

• `_range` = domain of integration

• `_expr` = integrand expression

Optional parameters:

• `_quad` = quadrature to use instead of the default one. Several ways are possible to pass the quadrature order for backward compatibility

 API Example Explanation Version `_quad=` `_quad=5` Pass the quadrature order as an integer, the default quadrature is used v0.105 `_quad=im()` `_quad=im(5)` Pass the default quadrature formula on a triangle to integrate exactly order 5 polynomials v0.105 `_quad=im(mesh,)` `_quad=im(mesh,5)` Pass the default quadrature formula on a mesh `mesh` to integrate exactly order 5 polynomials v0.105 `_quad=_Q()` `_quad=_Q<5>()` Pass the quadrature order at compile time to integrate exactly order 5 polynomials. It will be deprecated in a future release up to v0.105
 Starting from v0.105, quadratures are built at runtime and no more at compile time which means that quadrature orders can be adjusted dynamically, e.g from a command line option
• `_geomap` = type of geometric mapping to use, that is to say:

 Feel Parameter Description `GEOMAP_HO` High order approximation (same of the mesh) `GEOMAP_OPT` Optimal approximation: high order on boundary elements order 1 in the interior `GEOMAP_01` Order 1 approximation (same of the mesh)

### 1.2. Example

From `doc/manual/tutorial/dar.cpp`

``````form1( ... ) = integrate( _range = elements( mesh ),
_expr = f*id( v ) );``````

From `doc/manual/tutorial/myintegrals.cpp`

``````  // compute integral f on boundary
double intf_3 = integrate( _range = boundaryfaces( mesh ),
_expr = f );``````

From `doc/manual/advection/advection.cpp`

``````  // using default quadrature
form2( _test = Xh, _trial = Xh, _matrix = D ) +=
integrate( _range = internalfaces( mesh ),
_expr = ( averaget( trans( beta )*idt( u ) ) * jump( id( v ) ) )
+ penalisation*beta_abs*( trans( jumpt( trans( idt( u ) )) )
*jump( trans( id( v ) ) ) ),
_geomap = geomap );
// use deprecated _Q
form2( _test = Xh, _trial = Xh, _matrix = D ) +=
integrate( _range = internalfaces( mesh ),
_expr = ( averaget( trans( beta )*idt( u ) ) * jump( id( v ) ) )
+ penalisation*beta_abs*( trans( jumpt( trans( idt( u ) )) )
*jump( trans( id( v ) ) ) ),
_geomap = geomap );``````

From `doc/manual/laplacian/laplacian.cpp`

`````` auto l = form1( _test=Xh, _vector=F );
l = integrate( _range = elements( mesh ),
_expr=f*id( v ) ) +
integrate( _range = markedfaces( mesh, "Neumann" ),
_expr = nu*gradg*vf::N()*id( v ) );``````

## 2. Computing my first Integrals

This part explains how to integrate on a mesh with Feel++ (source `doc/manual/tutorial/myintegrals.cpp` ).

Let’s consider the domain $\Omega=[0,1]^d$ and associated meshes. Here, we want to integrate the following function

\begin{aligned} f(x,y,z) = x^2 + y^2 + z^2 \end{aligned}

on the whole domain $\Omega$ and on part of the boundary $\Omega$.

There is the appropriate code:

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

// create the mesh (specify the dimension of geometric entity)
auto mesh = unitHypercube<3>();

// our function to integrate
auto f = Px()*Px() + Py()*Py() + Pz()*Pz();

// compute integral of f (global contribution)
double intf_1 = integrate( _range = elements( mesh ),
_expr = f ).evaluate()( 0,0 );

// compute integral of f (local contribution)
double intf_2 = integrate( _range = elements( mesh ),
_expr = f ).evaluate(false)( 0,0 );

// compute integral f on boundary
double intf_3 = integrate( _range = boundaryfaces( mesh ),
_expr = f ).evaluate()( 0,0 );

std::cout << "int global ; local ; boundary" << std::endl
<< intf_1 << ";" << intf_2 << ";" << intf_3 << std::endl;
}``````

## 3. Mean value of a function

Let $f$ a bounded function on domain $\Omega$. You can evaluate the mean value of a function thanks to the `mean()` function :

$\bar{f}=\frac{1}{|\Omega|}\int_\Omega f=\frac{1}{\int_\Omega 1}\int_\Omega f$

### 3.1. Interface

``  mean( _range, _expr, _quad, _geomap );``

Required parameters:

• `_range` = domain of integration

• `_expr` = mesurable function

Optional parameters:

• `_quad` = quadrature to use.

• Default = `integer` corresponding to the quadrature order, see integrate.

• `_geomap` = type of geometric mapping.

• Default = `GEOMAP_OPT`

### 3.2. Example

Stokes example using `mean`
``````int main(int argc, char**argv )
{
Environment env( _argc=argc, _argv=argv,
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));

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

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

// element U=(u,p) in Vh
auto U = Vh->element();
auto u = U.element<0>();
auto p = U.element<1>();

// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),

a+= integrate(_range=elements(mesh),
_expr=-div(u)*idt(p)-divt(u)*id(p));

auto syms = symbols<2>();
auto u1 = parse( option(_name="functions.alpha").as<std::string>(), syms );
auto u2 = parse( option(_name="functions.beta").as<std::string>(), syms );
matrix u_exact = matrix(2,1);
u_exact = u1,u2;
auto p_exact = parse( option(_name="functions.gamma").as<std::string>(), syms );
auto f = -laplacian( u_exact, syms ) + grad( p_exact, syms ).transpose();
LOG(INFO) << "rhs : " << f;

// right hand side
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
_expr=trans(expr<2,1,5>( f, syms ))*id(u));
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
_expr=expr<2,1,5>(u_exact,syms));

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

double mean_p = mean(_range=elements(mesh),_expr=idv(p))(0,0);
double mean_p_exact = mean(_range=elements(mesh),_expr=expr(p_exact,syms))(0,0);
double l2error_u = normL2( _range=elements(mesh), _expr=idv(u)-expr<2,1,5>( u_exact, syms ) );
double l2error_p = normL2( _range=elements(mesh), _expr=idv(p)-mean_p-(expr( p_exact, syms )-mean_p_exact) );
LOG(INFO) << "L2 error norm u: " << l2error_u;
LOG(INFO) << "L2 error norm p: " << l2error_p;

// save results
auto e = exporter( _mesh=mesh );
e->save();
}``````

## 4. Norms

Let $f$ a bounded function on domain $\Omega$.

### 4.1. L2 norms

Let $f \in L^2(\Omega)$ you can evaluate the $L^2$ norm using the `normL2()` function:

$\parallel f\parallel_{L^2(\Omega)}=\sqrt{\int_\Omega |f|^2}$

#### 4.1.1. Interface

``  normL2( _range, _expr, _quad, _geomap );``

or squared norm:

``  normL2Squared( _range, _expr, _quad, _geomap );``

Required parameters:

• `_range` = domain of integration

• `_expr` = mesurable function

Optional parameters:

• `_quad` = quadrature to use.

• Default = `_Q<integer>()`

• `_geomap` = type of geometric mapping.

• Default = `GEOMAP_OPT`

#### 4.1.2. Example

From `doc/manual/laplacian/laplacian.cpp`

``````  double L2error =normL2( _range=elements( mesh ),
_expr=( idv( u )-g ) );``````

From `doc/manual/stokes/stokes.cpp`

Stokes example using `mean`
``````int main(int argc, char**argv )
{
Environment env( _argc=argc, _argv=argv,
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));

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

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

// element U=(u,p) in Vh
auto U = Vh->element();
auto u = U.element<0>();
auto p = U.element<1>();

// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),

a+= integrate(_range=elements(mesh),
_expr=-div(u)*idt(p)-divt(u)*id(p));

auto syms = symbols<2>();
auto u1 = parse( option(_name="functions.alpha").as<std::string>(), syms );
auto u2 = parse( option(_name="functions.beta").as<std::string>(), syms );
matrix u_exact = matrix(2,1);
u_exact = u1,u2;
auto p_exact = parse( option(_name="functions.gamma").as<std::string>(), syms );
auto f = -laplacian( u_exact, syms ) + grad( p_exact, syms ).transpose();
LOG(INFO) << "rhs : " << f;

// right hand side
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
_expr=trans(expr<2,1,5>( f, syms ))*id(u));
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
_expr=expr<2,1,5>(u_exact,syms));

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

double mean_p = mean(_range=elements(mesh),_expr=idv(p))(0,0);
double mean_p_exact = mean(_range=elements(mesh),_expr=expr(p_exact,syms))(0,0);
double l2error_u = normL2( _range=elements(mesh), _expr=idv(u)-expr<2,1,5>( u_exact, syms ) );
double l2error_p = normL2( _range=elements(mesh), _expr=idv(p)-mean_p-(expr( p_exact, syms )-mean_p_exact) );
LOG(INFO) << "L2 error norm u: " << l2error_u;
LOG(INFO) << "L2 error norm p: " << l2error_p;

// save results
auto e = exporter( _mesh=mesh );
e->save();
}``````

### 4.2. H^1 norm

In the same idea, you can evaluate the H1 norm or semi norm, for any function $f \in H^1(\Omega)$:

\begin{aligned} \parallel f \parallel_{H^1(\Omega)}&=\sqrt{\int_\Omega |f|^2+|\nabla f|^2}\\ &=\sqrt{\int_\Omega |f|^2+\nabla f * \nabla f^T}\\ |f|_{H^1(\Omega)}&=\sqrt{\int_\Omega |\nabla f|^2} \end{aligned}

where $*$ is the scalar product $\cdot$ when $f$ is a scalar field and the frobenius scalar product $:$ when $f$ is a vector field.

#### 4.2.1. Interface

``  normH1( _range, _expr, _grad_expr, _quad, _geomap );``

or semi norm:

``  normSemiH1( _range, _grad_expr, _quad, _geomap );``

Required parameters:

• `_range` = domain of integration

• `_expr` = mesurable function

• `_grad_expr` = gradient of function (Row vector!)

Optional parameters:

• `_quad` = quadrature to use.

• Default = `_Q<integer>()`

• `_geomap` = type of geometric mapping.

• Default = `GEOMAP_OPT`

normH1() returns a float containing the H^1 norm.

#### 4.2.2. Example

With expression:

``````  auto g = sin(2*pi*Px())*cos(2*pi*Py());
-2*pi*sin(2*pi*Px())*sin(2*pi*Py())*oneY();
// There gradg is a column vector!
// Use trans() to get a row vector
double normH1_g = normH1( _range=elements(mesh),
_expr=g,

With test or trial function `u`

``````  double errorH1 = normH1( _range=elements(mesh),
_expr=(u-g),

### 4.3. $L^\infty$ norm

You can evaluate the infinity norm using the normLinf() function:

$\parallel f \parallel_\infty=\sup_\Omega(|f|)$

#### 4.3.1. Interface

``  normLinf( _range, _expr, _pset, _geomap );``

Required parameters:

• `_range` = domain of integration

• `_expr` = mesurable function

• `_pset` = set of points (e.g. quadrature points)

Optional parameters:

• `_geomap` = type of geometric mapping.

• Default = `GEOMAP_OPT`

The `normLinf()` function returns not only the maximum of the function over a sampling of each element thanks to the `_pset` argument but also the coordinates of the point where the function is maximum. The returned data structure provides the following interface

• `value()`: return the maximum value

• `operator()()`: synonym to `value()`

• `arg()`: coordinates of the point where the function is maximum

#### 4.3.2. Example

``````  auto uMax = normLinf( _range=elements(mesh),
_expr=idv(u),
_pset=_Q<5>() );
std::cout << "maximum value : " << uMax.value() << std::endl
<<  "         arg : " << uMax.arg() << std::endl;``````