DUNE-FEM (unstable)
A fundamental task in any PDE solver is integration of a function over a domain given by a grid:
\[ I = \int_\Omega u(x)\,dx = \sum_{e\in E^0} \int_e u(x)\,dx \]
Here we demonstrate how this can be done for a user-defined function in global coordinates.
As always we need to set up a grid first. Here we use the Dune::YaspGrid class in four dimensions in order to demonstrate that even dimension larger than three can be used.
The Dune grid interface employs the class templates Dune::FieldVector and Dune::FieldMatrix to for small vectors and matrices of compile-time known size in many places, for example to represent positions and Jacobians of certain maps. Here is a small tutorial in the usage of these classes:
Next we define the function to integrate as a generic lambda function. We expect the argument x to be an instance of type Dune::FieldVector:
First we demonstrate how to compute an approximation to the integral by using the simple midpoint rule. This means traversing all elements, evaluate the function at the barycenter of the element, multiplying with the volume of the element and summing up:
A more accurate approximation of the integral for sufficiently smooth functions can be computed by quadrature rules. These are available in Dune for all the different element types and various integration orders. This simply requires an additional loop over the quadrature points within an element:
For integrating the divergence of a vector field \(f\) we might use the following formula:
\[ I = \int_\Omega \nabla\cdot f(x)\,dx = \int_{\partial \Omega} f(x)\cdot n(x) \,ds = \sum_{e\in E^0} \int_{\partial e\cap \partial \Omega} f(x)\cdot n(x) \,ds. \]
This is implemented by the following snippet illustrating the use of intersections:
Full example code: recipe-integration.cc