How to use block matrices and vectors?

When using the default BasisCreator provided by AMDiS, like LagrangeBasis or TaylorHoodBasis, the resulting indexing scheme is flat, meaning we need to use one (sparse) matrix and one vector to describe the linear system. Sometimes it is useful to have incorporate the blocking structure of a basis into the linear algebra.

Think of the stokes equation with velocity and pressure part. It could be useful to directly access the block of the matrix implementing the velocity Laplacian.

The blocking of the linear algebra matrix and vector is defined by the blocking structure of the global basis. In a blocked basis when accessing the index of a local function it returns not a integer, but a tuple of indices in referring to an entry in the hierarchy of vectors/matrices. See The interface for functions in the dune-functions module for some details about the indexing a index-merging strategies.

In order to use blocked matrices and vectors in AMDiS, you simply have to set the blocking strategy for the global basis. Since this is not done for the default BasisCreators, you need to create the basis directly using the factory functions from dune-functions:

#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>

#include <amdis/ProblemStat.hpp>

// the factory functions are defined in this namespace
using namespace Dune::Functions::BasisFactory;

// create a Taylor-Hood basis with blocking
auto gv = grid->leafGridView();
auto basis = makeBasis(gv,
  composite(
    power<2>(lagrange<2>(), blockedLexicographic()),
    lagrange<2>(), blockedLexicographic()
  ));

// Create a Problem
ProblemStat prob("prob", *grid, basis);
prob.initialize(INIT_ALL);

Now, the system matrix and system vector is not flat anymore, but blocked.

// reference to the internal vector structure
auto& vec = prob.solutionVector()->vector();

// we need integral constants to access vector components
using namespace Dune::Indices;

// the velocity and pressure block vector
auto& vec_vel = vec[_0];
auto& vec_pressure = vec[_1];

// the components of the velocity
auto& vec_vel_x = vec_vel[0];
auto& vec_vel_y = vec_vel[1];

Note: Currently the blocking is implemented for the ISTL backend only.