Skip to content

WIP: hierarchicVectorWrapper and hierarchicMatrixWrapper based on blocking structure

Simon Praetorius requested to merge simon.praetorius/dune-functions:feature/hierarchic_matrix_wrapper into master

This MR should provide the code for a discussion on this topic.

Summary

This MR adds a hierarchicVectorWrapper and hierarchicMatrixWrapper that provide MultiIndex access and resize functionality for hierarchically nested vectors and matrices.

Motivation

To use the blocking structure of the NodeFactories in linear-algebra data-structures, one has to provide a way to access the vector-coefficients and matrix-entries using a multi-index. This access requires a recursive call to access operators, like operator[.] for vectors and operator[.][.] (or operator(.,.)) for matrices (where the first one is a sequence of operator[] calls for row and columns indices.

The current implementation of a hierarchicVectorWrapper implements this recursion, using a common coefficient type as break-condition. This works fine in many cases, but fails if the coefficient is a wrapper type, like for std::vector<bool>. In the case of a hierarchicMatrixAccess, a similar strategy would be difficult to implement, due to the two-dimensional recursion. The coefficient type would be reached only if row and column recursion are finished. But it can not be checked separately.

This MR thus implements a different strategy. Here, the blocking-structure of the basis must be known. This blocking-structure is used as a break-condition for the recursions.

Implementation

Nodefactories can be equipped with the following IndexMergingStrategies:

  • Flat: FlatLexicographic, FlatInterleaved
  • Blocked: BlockedLexicographic
  • LeafBlocked: LeafBlockedInterleaved

So, it contains a blocking-structure and an indexing pattern. Only the blocking-structure is relevant for the multiindex-access. For ease of specialization on these block-structures, a light-weight representation is extracted (see functionspacebases/blocking.hh)

In the recursion of multiindex-access, the multiindex-component is increased while accessing the sub-blocking-structure. Thus, a break-condition is that the Blocking::Flat is reached (or, as a special case, Blocking::LeafBlocked).

Discussion

In order to provide a generic implementation, various hybrid-methods had to be implemented: hybridSize(), hybridNumRows(), and hybridNumCols(). These are hybrid in the sense, that is provides a size_t in case the vector (matrix) has dynamic index acces, and an index_constant otherwise. Thus, it can be used directly in the range() function of hybridutilities and in forEach() to have by default a (dynamic) loop and only if this is not possible, a static loop. I'm not sure whether this can be implemented in terms of other already existing functions in dune-common.

Different vector (matrix) types have different functions for resizing, e.g. a std::vector has a method resize(), while an MTL4 dense_vector uses change_dim() for resizing. It's even worth for matrices, i.e. some matrix-data-structure uses resize(.,.), another uses setSize(.,.) and the third one change_dim(.,.). There are two possible implementations for a generic interface:

  1. enforce one resize function. All data-structure that do not provide this, must be wrapped.
  2. provide a helper-class that does the resizing, templatized with the matrix-type. This can be specialized for another data-structure with another interface. (function overloading can not be used, i.e. because of the diffuse name-lookup rules the order of declaration and namespaces matter and nonetheless the wrong functions may be called. This is not the case for class specialization)

I followed the second choice, but it can be discussed whether possibility 1. should be considered instead.

Status

Not all tests are implemented yet. But for the hierarchicVectorWrapper one case fails that works with the old implementation:

static const std::size_t dim = 3;
using VelocityVector = std::vector<MultiTypeBlockVector<FieldVector<double,1>, double, FieldVector<double,1>>>;
using PressureVector = BlockVector<FieldVector<double,1>>;
using Vector = MultiTypeBlockVector<VelocityVector, PressureVector>;

because of the mixed types in the LeafBlocked velocity components. I think this is a very strange choice of data-structure, but maybe one could find a way to support this as well.

For the testing I have added a basis generator for several variants of a Taylor-Hood basis (different index-merging strategies). Maybe one could add more bases.

Merge request reports