Grid intersections

Intersections represent codimension-1 interfaces between neighboring elements or between an element and the domain boundary. They are central in finite-volume fluxes, DG methods, and boundary integrals.

Learning goals

After this page, you should be able to:

• iterate over intersections of each element,
• distinguish interior and boundary intersections,
• use intersection geometry and measures in practical loops.

Why intersections matter

Many PDE operators are defined on faces, not only on cell volumes. Typical examples are:

• boundary conditions,
• numerical fluxes across cell interfaces,
• jump/average terms in DG methods.

Iterating over intersections

Intersections are obtained per element:

for (const auto& e : elements(gv))
  for (const auto& is : intersections(gv, e))
  {
    if (is.boundary()) {
      // boundary face of e
    } else if (is.neighbor()) {
      // interior face between inside() and outside()
    }
  }

Useful methods:

• is.boundary(): true for domain boundary intersections,
• is.neighbor(): true if an outside element exists,
• is.inside() / is.outside(): adjacent elements,
• is.geometry(): geometry of the codim-1 intersection entity.

Intersection geometry

The geometry interface is analogous to element geometry:

const auto geo = is.geometry();
const auto center = geo.center();
const double measure = geo.volume();

In 2D, geo.volume() is edge length. In 3D, geo.volume() is face area.

Counting interior and boundary intersections

If you loop over all elements, interior intersections appear twice (once from each neighboring element). A common pattern is to keep only one side:

const auto& indexSet = gv.indexSet();
if (is.neighbor()) {
  auto in = indexSet.index(is.inside());
  auto out = indexSet.index(is.outside());
  if (in < out) {
    // count this interior intersection once
  }
}

Runnable example: examples/06_grid_intersections/intersections.cc

#include <array>
#include <iomanip>
#include <iostream>

#include <dune/common/fvector.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{4, 3};
  Dune::YaspGrid<dim> grid(L, N);
  const auto gv = grid.leafGridView();
  const auto& indexSet = gv.indexSet();

  std::size_t boundaryIntersections = 0;
  std::size_t interiorIntersections = 0;
  double boundaryMeasure = 0.0;
  double interiorMeasure = 0.0;

  for (const auto& e : elements(gv)) {
    for (const auto& is : intersections(gv, e)) {
      if (is.boundary()) {
        ++boundaryIntersections;
        boundaryMeasure += is.geometry().volume();
      } else if (is.neighbor()) {
        const auto insideIndex = indexSet.index(is.inside());
        const auto outsideIndex = indexSet.index(is.outside());
        if (insideIndex < outsideIndex) {
          ++interiorIntersections;
          interiorMeasure += is.geometry().volume();
        }
      }
    }
  }

  std::cout << "elements: " << gv.size(0) << "\n";
  std::cout << "boundary intersections: " << boundaryIntersections << "\n";
  std::cout << "interior intersections: " << interiorIntersections << "\n";
  std::cout << std::fixed << std::setprecision(6);
  std::cout << "boundary measure: " << boundaryMeasure << "\n";
  std::cout << "interior measure: " << interiorMeasure << "\n";

  return 0;
}

Build and run:

source _env/activate.sh
cd examples/06_grid_intersections
dunecontrol --current all
./build/debug/example06

Output:

elements: 12
boundary intersections: 14
interior intersections: 17
boundary measure: 4.000000
interior measure: 5.000000

Summary and next steps

You now have the face-level traversal pattern used for fluxes and boundary terms. Continue with Attaching data to grid entities, then Grid hierarchy and refinement.