How to use block matrices and vectors?
When using the default
BasisCreator provided by AMDiS, like
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
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
#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; auto& vec_vel_y = vec_vel;
Note: Currently the blocking is implemented for the ISTL backend only.