Feel++ 0.91.0
|
We would like to compute some integrals on a domain and parts of the domain, i.e. subregions and (parts of) boundary. In this tutorial the domain can be either
or
with and
.
These domains are plotted in this Section.
Once we have defined the computational mesh, we would like to compute the following integrals:
We now turn to the implementation. To compute the first integral, the code reads
//# marker1 # double local_domain_area = integrate( elements(mesh), constant(1.0)).evaluate()(0,0); //# endmarker1 #
elements(mesh)
returns a pair of iterators over the elements owned by the current processor. Now that we have computed the integral of 1 over the region of current processor (ie the area of the domain owned by the processor), we want to compute the area of
. To do that we collect the integrals on all processors using a
reduce
MPI operation and sum all contributions. We have used here the Boost.MPI library that provides an extremely powerful wrapper around the MPI library. The code reads
//# marker2 # double global_domain_area=local_domain_area; if ( this->comm().size() > 1 ) mpi::all_reduce( this->comm(), local_domain_area, global_domain_area, std::plus<double>() ); //# endmarker2 #
Finally, we print to the log file the result of the local and global integral calculation.
//# marker3 # Log() << "int_Omega 1 = " << global_domain_area << "[ " << local_domain_area << " ]\n"; //# endmarker3 #
Another calculation is for example to compute the perimeter of the domain
//# marker4 # double local_boundary_length = integrate( boundaryfaces(mesh), constant(1.0)).evaluate()(0,0); double global_boundary_length = local_boundary_length; if ( this->comm().size() > 1 ) mpi::all_reduce( this->comm(), local_boundary_length, global_boundary_length, std::plus<double>() ); Log() << "int_BoundaryOmega (1)= " << global_boundary_length << "[ " << local_boundary_length << " ]\n"; //# endmarker4 #
The main difference with the domain area computation resides in the elements with are iterating on: here we are iterating on the boundary faces of the domain to compute the integral using boundaryfaces(mesh)
to provide the pairs of iterators.
Now say that we want to compute the third integral. The Finite Element Embedded Language (FEEL++) language provides the keyword Px()
and Py()
to denote the and
coordinates like in equation above. The code reads then
//# marker5 # double local_intf = integrate( elements(mesh), Px()*Px() + Py()*Py() + Pz()*Pz() ).evaluate()(0,0); //# endmarker5 #
Feel computes automatically the quadrature associated the expression under the integral. If the expression is polynomial then the quadrature is exact. If the expression is not polynomial, then each non-polynomial term in the expression is considered a polynomial of degree 2 by default.
Here is how the 4-th integral can be computed letting Feel decide about the quadrature
//# marker6 # double local_intsin = integrate( elements(mesh), sin( Px()*Px() + Py()*Py() + Pz()*Pz() ) ).evaluate()(0,0); //# endmarker6 #
An alternative is to set yourself the quadrature by passing a new argument to the integrate
function: _Q<Order>()
where Order
is the maximum polynomial order that the quadrature should integrate exactely.
//# marker7 # double local_intsin2 = integrate( elements(mesh), _Q<2>(), sin( Px()*Px() + Py()*Py() + Pz()*Pz() ) ).evaluate()(0,0); //# endmarker7 #
Let's run now the tutrial example feel_doc_myintegrals
: First we generate the simplex shape domain in all 3 dimensions
feel_doc_myintegrals --shape="simplex"--nochdir
![]() Line | ![]() Triangle | ![]() Tetrahedron |
The output gives
We do the same now with the hypercube shape domain in all 3 dimensions
feel_doc_myintegrals --shape="hypercube"--nochdir
![]() Line | ![]() Unit Square | ![]() Unit Cube |
The output gives