Skip to content

Integrating functions on a grid

Integration is one of the most common tasks in PDE codes. This page shows two practical approaches based on dune-grid and dune-geometry.

Learning goals

After this page, you should be able to:

  • implement element-wise integration over a grid,
  • compare midpoint and quadrature-based integration,
  • use intersections for boundary contributions.

Problem

Given a function \(u\colon\Omega\to\mathbb{R}\), we want to compute the integral \(\int_\Omega u(x)\,dx\) over a grid domain \(\Omega\).

1. Midpoint rule

Use this as a simple and fast baseline integration method.

double integral = 0.0;
auto gv = grid.leafGridView();

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

for (const auto& e : elements(gv))
{
  const auto geo = e.geometry();
  integral += u(geo.center()) * geo.volume();
}

This is easy and cheap, but only low-order accurate.

2. Quadrature rule

Use this for the standard finite-element integration workflow.

#include <dune/geometry/quadraturerules.hh>

double integral2 = 0.0;
constexpr int dim = Grid::dimension;
using QR = Dune::QuadratureRules<typename Grid::ctype, dim>;

for (const auto& e : elements(gv))
{
  const auto geo = e.geometry();
  const auto quadrature = QR::rule(geo.type(), 5);  // order 5
  for (const auto& qp : quadrature)
    integral2 += u(geo.global(qp.position()))
               * geo.integrationElement(qp.position())
               * qp.weight();
}

This is the typical pattern used in element assembly.

3. Boundary fluxes via intersections

For boundary integrals, iterate over intersections and test boundary sides:

double flux = 0.0;
auto f = [](const auto& x) { return x; };

for (const auto& e : elements(gv))
  for (const auto& is : intersections(gv, e))
    if (is.boundary())
      flux += f(is.geometry().center()) * is.centerUnitOuterNormal() * is.geometry().volume();

For a focused explanation of this interface, see Grid intersections.

Runnable example: examples/03_grid_integration/integration.cc 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>

#include <dune/common/fvector.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/grid/yaspgrid.hh>

int main()
{
  constexpr int dim = 2;
  Dune::FieldVector<double, dim> L(1.0);
  std::array<int, dim> N{16, 16};
  Dune::YaspGrid<dim> grid(L, N);
  const auto gv = grid.leafGridView();

  auto u = [](const auto& x) {
    return x[0] * x[0] + x[1];
  };

  double midpointIntegral = 0.0;
  for (const auto& e : elements(gv)) {
    const auto geo = e.geometry();
    midpointIntegral += u(geo.center()) * geo.volume();
  }

  double quadratureIntegral = 0.0;
  using QR = Dune::QuadratureRules<double, dim>;
  for (const auto& e : elements(gv)) {
    const auto geo = e.geometry();
    const auto quadrature = QR::rule(geo.type(), 4);
    for (const auto& qp : quadrature) {
      quadratureIntegral += u(geo.global(qp.position()))
                            * geo.integrationElement(qp.position())
                            * qp.weight();
    }
  }

  const double exact = 5.0 / 6.0;
  std::cout << std::fixed << std::setprecision(10);
  std::cout << "elements: " << gv.size(0) << "\n";
  std::cout << "midpoint integral  = " << midpointIntegral << "\n";
  std::cout << "quadrature integral= " << quadratureIntegral << "\n";
  std::cout << "exact integral     = " << exact << "\n";
  std::cout << "midpoint error     = " << std::abs(midpointIntegral - exact) << "\n";
  std::cout << "quadrature error   = " << std::abs(quadratureIntegral - exact) << "\n";

  return 0;
}

Build and run:

source _env/activate.sh
cd examples/03_grid_integration
dunecontrol --current all
./build/debug/example03

Output:

elements: 256
midpoint integral  = 0.8330078125
quadrature integral= 0.8333333333
exact integral     = 0.8333333333
midpoint error     = 0.0003255208
quadrature error   = 0.0000000000

Source inspiration

This page is based on the DUNE grid recipe: _material/dune-grid-recipes/recipe-integration.cc.

Summary and next steps

You now have the core integration loop patterns used in finite element assembly. Next, continue with Grid I/O and Grid info utilities.