Integrate

Feel++ provides the integrate() function to define integral expressions which can be used to compute integrals and to build linear or bilinear forms.

1. Interface

integrate( _range, _expr, _quad, _geomap );

The parameters are Boost-style named parameters; order is not important when names are used. Two parameters are required and two are optional.

1.1. Required parameters

  • _range — domain of integration (mesh iterators)

  • _expr — integrand expression

1.2. Optional parameters

  • _quad — quadrature to use. Examples:

API

Example

Explanation

Version

_quad=<n>

_quad=5

Pass the quadrature order as an integer

v0.105

_quad=im<Convex>(<n>)

_quad=im<Triangle>(5)

Use the default quadrature on triangles for polynomials up to order 5

v0.105

_quad=im(mesh,<n>)

_quad=im(mesh,5)

Use the mesh’s default quadrature of order 5

v0.105

_quad=_Q<n>()

_quad=_Q<5>()

Compile-time quadrature spec (deprecated as of v0.105)

up to v0.105
From v0.105 quadratures are built at runtime so quadrature orders can be adjusted dynamically.

  • _geomap — geometric mapping type. Options: GEOMAP_HO, GEOMAP_OPT, GEOMAP_01.

2. Examples

From doc/manual/tutorial/dar.cpp:

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

Compute an integral on the boundary:

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

Using a custom quadrature and geomap in a face-based integral:

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 );

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) );

3. Quick tutorial: compute a simple integral

Consider the domain \Omega = [0,1]^d and the integrand f(x,y,z) = x2+y2+z^2.

auto mesh = unitHypercube<3>();
auto f = Px()*Px() + Py()*Py() + Pz()*Pz();

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

// local (per-element) contributions
double intf_2 = integrate( _range = elements(mesh), _expr = f ).evaluate(false)(0,0);

// boundary integral
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;