Attaching data to grid entities

Most numerical codes need arrays of values attached to grid entities (for example per-cell coefficients, per-face fluxes, or per-vertex data). This page introduces index sets, id sets, and mappers for this task.

Learning goals

After this page, you should be able to:

• use indexSet to store data in contiguous vectors,
• distinguish indexSet from idSet,
• use MultipleCodimMultipleGeomTypeMapper for mixed geometry types and multiple codimensions.

Index set vs id set

For a grid view gv:

const auto& indexSet = gv.indexSet();
const auto& idSet = grid.globalIdSet();

• indexSet gives consecutive, zero-based indices for entities in a view,
• idSet gives globally unique ids, intended to stay persistent under grid changes.

Important detail: indexSet numbering is defined per codimension and geometry-type class. In mixed 2D grids, triangles and quads are therefore numbered independently. If you need one contiguous numbering across both types, use a mapper.

Typical size queries with indexSet:

indexSet.size(0);    // elements
indexSet.size(1);    // codim-1 entities
indexSet.size(dim);  // vertices

Data attached with indexSet

For edge-based data in 2D (codim 1):

std::vector<double> edgeLengths(indexSet.size(1), 0.0);
for (const auto& edge : edges(gv))
  edgeLengths[indexSet.index(edge)] = edge.geometry().volume();

This is the standard fast storage pattern for many assembly/data loops.

Data attached with a mapper

MultipleCodimMultipleGeomTypeMapper creates contiguous indices for selected subsets. For element-based data:

using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
Mapper elementMapper(gv, Dune::mcmgElementLayout());

std::vector<double> cellValue(elementMapper.size(), 0.0);
for (const auto& e : elements(gv))
  cellValue[elementMapper.index(e)] = e.geometry().center().two_norm();

This pattern is useful if the subset is more specific than plain codim sizes.

One contiguous container over multiple codimensions

Mappers are also useful if one container should hold data for several codimensions, for example elements and edges together:

const Dune::MCMGLayout elementAndEdgeLayout = [](Dune::GeometryType gt, int griddim) {
  const auto gridDim = static_cast<unsigned int>(griddim);
  return (gt.dim() == gridDim) || (gt.dim() + 1 == gridDim);
};

Mapper mixedMapper(gv, elementAndEdgeLayout);
std::vector<double> mixedData(mixedMapper.size(), 0.0);

Here mixedMapper.index(entity) maps all selected entities into a single contiguous index range, independent of codimension and geometry type.

Runnable example: examples/07_grid_data/grid_data.cc

#include <array>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <set>
#include <vector>

#include <dune/common/fvector.hh>
#include <dune/grid/common/mcmgmapper.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);
  auto gv = grid.leafGridView();

  const auto& indexSet = gv.indexSet();
  const auto& idSet = grid.globalIdSet();

  // Example 1: data attached to edges using indexSet indices (codim 1 in 2D).
  std::vector<double> edgeLengths(indexSet.size(1), 0.0);
  for (const auto& edge : edges(gv))
    edgeLengths[indexSet.index(edge)] = edge.geometry().volume();
  const double totalEdgeLength =
      std::accumulate(edgeLengths.begin(), edgeLengths.end(), 0.0);

  // Example 2: data attached to elements using a mapper.
  using GridView = decltype(gv);
  using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
  Mapper elementMapper(gv, Dune::mcmgElementLayout());
  std::vector<double> cellValue(elementMapper.size(), 0.0);
  for (const auto& e : elements(gv))
    cellValue[elementMapper.index(e)] = e.geometry().center().two_norm();
  const double cellValueSum =
      std::accumulate(cellValue.begin(), cellValue.end(), 0.0);

  // Example 3: one contiguous container for codim-0 and codim-1 entities.
  const Dune::MCMGLayout elementAndEdgeLayout = [](Dune::GeometryType gt, int griddim) {
    const auto gridDim = static_cast<unsigned int>(griddim);
    return (gt.dim() == gridDim) || (gt.dim() + 1 == gridDim);
  };
  Mapper mixedMapper(gv, elementAndEdgeLayout);
  std::vector<double> mixedData(mixedMapper.size(), 0.0);
  for (const auto& e : elements(gv))
    mixedData[mixedMapper.index(e)] = 1.0;
  for (const auto& edge : edges(gv))
    mixedData[mixedMapper.index(edge)] = 2.0;
  const double mixedDataSum =
      std::accumulate(mixedData.begin(), mixedData.end(), 0.0);

  // Example 4: global ids are unique and persistent under grid modifications.
  using VertexId = typename std::remove_cvref_t<decltype(idSet)>::IdType;
  std::set<VertexId> vertexIds;
  for (const auto& v : vertices(gv))
    vertexIds.insert(idSet.id(v));

  std::cout << "elements: " << indexSet.size(0) << "\n";
  std::cout << "edges: " << indexSet.size(1) << "\n";
  std::cout << "vertices: " << indexSet.size(dim) << "\n";
  std::cout << "unique vertex ids: " << vertexIds.size() << "\n";
  std::cout << "mixed mapper size (elements + edges): " << mixedMapper.size() << "\n";
  std::cout << std::fixed << std::setprecision(6);
  std::cout << "sum of edge lengths: " << totalEdgeLength << "\n";
  std::cout << "sum of mapped cell values: " << cellValueSum << "\n";
  std::cout << "sum of mixed data: " << mixedDataSum << "\n";

  return 0;
}

Build and run:

source _env/activate.sh
cd examples/07_grid_data
dunecontrol --current all
./build/debug/example07

Output:

elements: 12
edges: 31
vertices: 20
unique vertex ids: 20
mixed mapper size (elements + edges): 43
sum of edge lengths: 9.000000
sum of mapped cell values: 9.111313
sum of mixed data: 74.000000

Summary and next steps

You now have the core data-attachment tools used in practical DUNE codes. Continue with Grid hierarchy and refinement, then Integrating functions on a grid.