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=<n>

_quad=5

Pass the quadrature order as an integer, the default quadrature is used

v0.105

_quad=im<Convex>(<n>)

_quad=im<Triangle>(5)

Pass the default quadrature formula on a triangle to integrate exactly order 5 polynomials

v0.105

_quad=im(mesh,<n>)

_quad=im(mesh,5)

Pass the default quadrature formula on a mesh mesh to integrate exactly order 5 polynomials

v0.105

_quad=_Q<n>()

_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 ),
               _quad = 2*Order,
               _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 ),
               _quad = _Q<2*Order>(),
               _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(),
                     _about=about( _name="myintegrals" ,
                                   _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,
                     _about=about(_name="mystokes",
                                  _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),
                  _expr=trace(gradt(u)*trans(grad(u))) );

    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->add( "u", u );
    e->add( "p", p );
    e->save();
}

4. Norms

Let \(f\) a bounded function on domain \(\Omega\).

4.1. L² 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,
                     _about=about(_name="mystokes",
                                  _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),
                  _expr=trace(gradt(u)*trans(grad(u))) );

    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->add( "u", u );
    e->add( "p", p );
    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());
  auto gradg = 2*pi*cos(2* pi*Px())*cos(2*pi*Py())*oneX()
            -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,
                   _grad_expr=trans(gradg) );

With test or trial function u

  double errorH1 = normH1( _range=elements(mesh),
                    _expr=(u-g),
                  _grad_expr=(gradv(u)-trans(gradg)) );

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;

Integration

1. Integrals

1.1. Interface

1.2. Example

2. Computing my first Integrals

3. Mean value of a function

3.1. Interface

3.2. Example

4. Norms

4.1. L2 norms

4.1.1. Interface

4.1.2. Example

4.2. H^1 norm

4.2.1. Interface

4.2.2. Example

4.3. \(L^\infty\) norm

4.3.1. Interface

4.3.2. Example

4.1. L² norms