The next step is to compute integrals over the mesh. The source code is available in myintegrals.cpp
.
Step by step explanations
We start by loading a Mesh in 2D
auto mesh = loadMesh( _mesh=new Mesh<Simplex<2>> );
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()
.
auto g = expr( soption(_name="functions.g") );
then We compute two integrals over the domain and its boundary respectively
- \(\int_{\partial \Omega} g\)
and we print the results to the screen.
auto intf_1 = integrate( _range = elements( mesh ),
_expr = g ).evaluate();
auto intf_2 = integrate( _range = boundaryfaces( mesh ),
_expr = g ).evaluate();
auto grad_g = grad<2>(g);
auto intgrad_f = integrate( _range = elements( mesh ),
_expr = grad_g ).evaluate();
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;
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.
shell> ./feelpp_doc_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
shell> ./feelpp_doc_myintegrals --functions.g=x:x
int_Omega x = 0.5
int_{boundary of Omega} x = 2
Recall that
: 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
shell> ./feelpp_doc_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_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
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
Complete code
The complete code reads as follows
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feelvf/integrate.hpp>
#include <feel/feelvf/ginac.hpp>
#include <feel/feelfilters/exporter.hpp>
int
main( int argc, char** argv )
{
Environment env( _argc=argc, _argv=argv,
_about=about( _name="myintegrals" ,
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org" ) );
auto mesh = loadMesh( _mesh=new Mesh<Simplex<2>> );
auto g = expr( soption(_name="functions.g") );
auto intf_1 = integrate( _range = elements( mesh ),
_expr = g ).evaluate();
auto intf_2 = integrate( _range = boundaryfaces( mesh ),
_expr = g ).evaluate();
auto grad_g = grad<2>(g);
auto intgrad_f = integrate( _range = elements( mesh ),
_expr = grad_g ).evaluate();
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;
}
to compile just type
make feelpp_doc_myintegrals