[dune-functions] [Dune-devel] hierarchicVectorWrapper and hierarchicMatrixWrapper

[Carsten Gräser]

Am 21.02.2018 um 13:42 schrieb Carsten Gräser:
> Hi Simon,
> we have a dedicated dune-functions mailing list
> 
>   dune-functions at lists.dune-project.org
> 
> and this discussion should probably go there.
> Hence I'll reply to this list regarding your
> proposal.
OK, here's my actual reply on the topic:

WARNING:
@joe-average-user: The following discussion is
quite technical and maybe not interesting to you.

@developers: Simon has a point here and I'm
also not happy with the current situation.

> Am 21.02.2018 um 13:17 schrieb Simon Praetorius:
>> Hi folks,
>>
>> during my development of a discretization module based directly on
>> dune-functions, I've implemented wrappers for hierarchically nested
>> vectors and matrices. The former was already available, but I don't see
>> how it can be extended from vectors to matrices. Thus, I've made a new
It can be done for matrices in (almost) exactly the same way.
I did it for dune-fufem. However, there's one issue when using
non-uniform multi-indices to access a matrix:

For non-diagonal entries you will query entries where the
row- and column-multi-index have a different length. In
this case some indices conceptually refer to a single-row
matrix that does not require a column index (or the other
way around). Since our matrix implementations do not support
this, I solved this in dune-fufem by providing wrappers that
essentially tag a matrix in the nested matrix structure to
be single-row or single-column. Then the recursive multi-index
resolver (working similar to the vector case) can smuggle
in a [0] for the missing index.

This is not in dune-functions, because some considered this
a feature of a discretization module that should thus not
be part of a dune-functions. I think that this is debatable
though.

>> approach and came up with the idea of coupling the hierarchically nested
>> vector/matrix with the corresponding nesting structure. When declaring
>> the vector/matrix type this must already be known and is naturally
>> related to the index-merging strategy of the GlobalBasis.
>>
>> This coupling of data-structure and nesting-structure works quite fine
>> and (nearly) all the former tests for hierarchicVectorWrapper can be
>> passed. Additionally, some more data-structures are supported, e.g. for
>> a taylor-hood basis
>>
>> - TupleVector<std::array<std::vector<bool>,dim>, std::vector<bool>>
>> - TupleVector<BitFieldVector<dim>, std::vector<bool>>
>>
>> usable as a bit-field in interpolation methods, for example.
>>
>> With "nesting-structure" I mean a light-weight representation of the
>> index-merging strategies, e.g. for the examples above:
>>
>> - Blocked<Blocked<Flat>, Flat>
>> - Blocked<LeafBlocked<dim>, Flat>
>>
>> where `Blocked` is related to `BlockedLexicographic`, `Flat` is related
>> to (`FlatLexicographic`, `FlatInterleaved`) and `LeafBlocked<dim>` is
>> related to `LeafBlockedInterleaved`.
>>
>> This development is now ported back from our discretization module to
>> dune-functions and is uploaded as an experimental branch, see MR !90.
>> Currently, a lot of utility code is necessary. But maybe this can be
>> eliminated or replaced with other code.
>>
>> What do you think? Is this a way to go or should we look for a another
>> way to write wrappers purely on the data-structure without knowledge of
>> the nesting-structure?
This looks interesting. Since one cannot terminate the
recursion depending on the dynamic length of the multi-index
I chose to terminate depending on the data structure.
While dune-functions is intended to work with scalar
coefficients (vector-valued bases can naturally be
implemented as power-*), prescribing the scalar
coefficient type is nasty. On the other hand I do not
see any possibility to avoid a common coefficient type
if you have dynamic multi-indices.

What I like about the approach you outlined, is that
it enriches the information about the index tree structure.
But hard-coding the Blocked/Flat information is unfortunate,
because dune-functions does not require that a basis
is implemented as DefaultGlobalBasis<PreBasis> (see
my recent mail on the renaming). Implementing your
own basis without this mechanism is totally fine.

I'd rather like to view this as some abstract
information about the index set additionally to
size(prefix). For DefaultLocalIndexSet it could
be implemented using the knowledge about the
IndexMergingStrategies. My rough idea about this
was to have something like isStatic(prefix)
and staticSize(prefix). I have the impression
that your proposal essentially does this, but I'll
have to look into the details.

BTW:
There's an alternative approach that I had in mind.
The whole problem arises from having dynamically
sized global multi-indices which are unavoidable
if we want to map our _flat_ local indices to global
multi-indices. If we additionally provide a method
index(treePath,index) with a static treePath to
a leaf and an index of the basis function only within
this leaf, then this could make use of the static
information and return a hybrid (static-dynamic)
version of the multi-index. Resolving this would
be a trivial recursion. But it would also require
to implement the loop in th local->global transformation
in user code to a tree traversal.

Anyway, thanks for your effort and for sharing your code.

Best,
Carsten

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 455 bytes
Desc: OpenPGP digital signature
URL: <https://lists.dune-project.org/pipermail/dune-functions/attachments/20180221/d6dc4687/attachment.sig>