DUNE-FEM (unstable)

Integrate a function on a grid

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.

const int dim = 4;
using Grid = Dune::YaspGrid<dim>;
Dune::FieldVector<double,dim> len; for (auto& l:len) l=1.0;
std::array<int,dim> cells; for (auto& c : cells) c=5;
Grid grid(len,cells);
vector space out of a tensor product of fields.
Definition: fvector.hh:91
[ provides Dune::Grid ]
Definition: yaspgrid.hh:166

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:

Dune::FieldVector<double,4> x({1,2,3,4}); // make a vector
auto y(x); // copy constructor
y *= 1.0/3.0; // scaling
[[maybe_unused]] auto s = x*y; // scalar product
[[maybe_unused]] auto norm = x.two_norm(); // Euclidean norm
Dune::FieldMatrix<double,4,4> A({{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}); // make a matrix
A.mv(x,y); // matvec: y = Ax
A.usmv(0.5,x,y); // axpy: y += 0.5*Ax
A dense n x m matrix.
Definition: fmatrix.hh:117

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:

auto u = [](const auto& x){return std::exp(x.two_norm());};

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:

double integral=0.0;
auto gv = grid.leafGridView(); // extract the grid view
for (const auto& e : elements(gv))
integral += u(e.geometry().center())*e.geometry().volume();
std::cout << "integral = " << integral << std::endl;
IteratorRange<... > elements(const GV &gv)
Iterates over all elements / cells (entities with codimension 0) of a GridView.

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:

double integral2 = 0.0;
for (const auto& e : elements(gv))
{
auto geo = e.geometry();
auto quadrature = QR::rule(geo.type(),5);
for (const auto& qp : quadrature)
integral2 += u(geo.global(qp.position()))
*geo.integrationElement(qp.position())*qp.weight();
}
std::cout << "integral2 = " << integral2 << std::endl;
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:260

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:

auto f = [](const auto& x){return x;};
double divergence=0.0;
for (const auto& i : elements(gv)) {
for (const auto& I : intersections(gv,i))
if (!I.neighbor())
{
auto geoI = I.geometry();
divergence += f(geoI.center())*I.centerUnitOuterNormal()*geoI.volume();
}
}
std::cout << "divergence = " << divergence << std::endl;
IteratorRange<... > intersections(const GV &gv, const Entity &e)
Iterates over all Intersections of an Entity with respect to the given GridView.

Full example code: recipe-integration.cc

Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 20, 23:30, 2024)