Step 1: Grid and linear algebra setup

This step creates the computational grid, global DoF numbering, and sparse matrix layout.

Learning goals

After this page, you should be able to:

• create a unit-square YaspGrid,
• define one global DoF per vertex using an index set,
• initialize BCRSMatrix and BlockVector.

Grid setup

For this initial example, we will use a structured grid on \((0,1)^2\) build from \(32\times 32\) square elements:

using Grid = Dune::YaspGrid<2>;
Grid grid({1.0, 1.0}, {32, 32});
auto gv = grid.leafGridView();

For a different domain shape, number of elements, element types, or even local refinement, we refer to grid introduction in Grid Basics.

Vertex DoF numbering

First-order conforming elements have one DoF attached to each vertex. Thus, we get as total number of degrees of freedom the number of vertices in the grid:

const auto numDofs = gv.indexSet().size(2);

In each element of the grid, the local corner vertices are associated to a global vertex index. This can be retrieved from the grid's index set, by collecting the vertex sub-indices for each element. For performance reasons, it is good to collect this index mapping of local to global vertex indices into a container,

template<class GridView, class Element>
std::vector<std::size_t> elementVertexDofs(const GridView& gv, const Element& e)
{
  constexpr int dim = GridView::dimension;
  const auto& indexSet = gv.indexSet();

  const int nCorners = e.geometry().corners();
  std::vector<std::size_t> dofs(nCorners);
  for (int localVertex = 0; localVertex < nCorners; ++localVertex)
    dofs[localVertex] = indexSet.subIndex(e, localVertex, dim);

  return dofs;
}

What this does:

• e.geometry().corners() gives the number of vertices of the current element (works in any dimension and for different geometry types).
• indexSet.subIndex(e, localVertex, dim) maps one local vertex of this element to the corresponding global index of the codim-dim entity, i.e. a vertex.
• The returned vector is the local-to-global map used when scattering local element contributions into global matrix/vector entries.

Note

If the number of corners or the number of DoFs per element is a statically known constant, one could store this mapping in a std::array, instead of a dynamic std::vector. This requires to know the element type of the grid and its dimension.

For example, a 2d YaspGrid allows only square elements, and a first order Lagrange local finite-element has one degree of freedom per vertex. This makes \(4=2^{\text{dim}}\) DoFs per Element.

constexpr int numElementDofs = 1<<dim;
std::array<std::size_t,numElementDofs> dofs;
for (int localVertex = 0; localVertex < numElementDofs; ++localVertex)
  dofs[localVertex] = indexSet.subIndex(e, localVertex, dim);

Sparse matrix pattern

Before assembly, we collect all nonzero couplings with MatrixIndexSet:

Dune::MatrixIndexSet pattern(numDofs, numDofs);

for (const auto& e : Dune::elements(gv)) {
  const auto dofs = elementVertexDofs(gv, e);
  for (std::size_t i = 0; i < dofs.size(); ++i)
    for (std::size_t j = 0; j < dofs.size(); ++j)
      pattern.add(dofs[i], dofs[j]);
}

Dune::BCRSMatrix<double> A;
pattern.exportIdx(A); // allocate row sizes + column indices from pattern
A = 0.0;

This pattern contains all row/column pairs that can receive element contributions. For a vertex-based first-order method, all local vertex DoFs on one element couple with each other, so we insert the full local nCorners x nCorners block.

The rhs will be assembled into a vector with numDofs entries,

Dune::BlockVector<double> b(numDofs);
b = 0.0;

Summary and next steps

Grid and linear algebra objects are ready for assembly. Continue with Step 2: Assembly and boundary conditions.