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

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

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, evalute 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;

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;
  using QR = Dune::QuadratureRules<Grid::ctype,dim>;
  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;

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;

Full example code: recipe-integration.cc