Misuse of [MatrixType]::mv breaks all transformed global bases

Currently all transformed bases use [MatrixType]::mv in a non-supported way. This breaks compilation of such bases:

• Hermite-type bases like MorleyBasis use FieldMatrix::mv in TransformedLocalBasis to transform to the element specific basis functions.
• Piola-transformed bases use jacobianInverseTransposed.mv(...) to transform the range space.

In both cases non-supported vector types like std::vector and std::array are passed to mv. This is not supported by the used matrix classes which breaks compilation. The reason that this was not discovered earlier is that the mv implementation normally works fine for these container types. But the new pipeline uses DUNE_FMatrix_WITH_CHECKING to activate runtime size checks which query the vector size using x.N() instead of x.size().

Notice that fixing this is not straight forward:

• If passing non-dune vectors to mv should be supported, then this is a bug in dune-common which should use x.size() instead of x.N().
• If mv requires that the passed vectors support the extended interface, then this should be fixed in dune-functions. Here we again have several options:
  • If we want to support passing raw vector types to interpolate as we currently do, then we either have to either copy the vectors before passing it to mv in the Piola-transformation or manually implement the mv operation locally. The latter should be avoided, because we don't have the matrix internals under control.
  • If we require that all vector types should support the extended interface (also for interpolate and function ranges), then this will break a lot of downstream code.

I don't see what's the problem in interpolate(). While we pass a container, this is not used in the local interpolation. Instead we create our own std::vector and pass it to the local finite-element. Did I miss another spot?

E.g. for Raviart-Thomas, interpolate has to wrap the passed function in order to apply the inverse contra-variant Piola transoform. The latter is implemented in terms of mtv. Here the argument vector type results from the range type of the user-provided function.

I have reported this some time ago, but then forgot during the merge about the issue. Recently I have provided a possible fix: !499. But this requires some utility in dune-common first. Essential, it allows to transform "any" (contiguous) vector-like type into the DenseVector interface. Then it can easily be passed to .mv and others.

core/dune-common!1354 is only a partial fix here, because we cannot expect that containers used in interpolate use contiguous storage. Before reading your comment I quickly implemented a DenseVectorView (similar to ScalarVectorView) that derives from DenseVector but uses a reference to a random access container. This can serve as a quick fix here.

But this is what I had already implemented some time ago and put there as a MR. A DenseVectorView (just called DenseVectorSpan because it only needs the underlying data handle - this could even be an iterator (iterators would need a slightly different interface)), that derives from DenseVector and has deduction guides for several std data structures and some Dune data structures to make life easier. This could always be used whenever we need to fulfill the DenseVector interface. Since it works on the raw data, there is no performance overhead. It even provides the proper static or dynamic size information automatically, thanks to the Std::[md]span base class.

I was not aware that std::mdpan can handle non-contiguous storage and proxy-containers.

Everything is handled by an accessor-type. This needs to be capable of return a reference to the ith element relative to an offset on a data handle. If the data handle is a random-access iterator it, it could easily return it[i]. mdspan is especially designed for strided laouts. There, the underlying data is contiguous, but is accessed only with strides.

But sure, if it is just about the .N(), a simpler wrapper that just adds this method to the interface of another container is simpler. But it is another of these artificial wrapper types (like also ScalarVectorView) one should get rid of eventually and design a proper interface.

Another obstacle is that DenseVectorSpan is not generic but needs dedicated deduction guides for all supported types.

one should get rid of eventually and design a proper interface.

This would require also require some automatic mechanism to deduce a suitable accessor and size for a generic container.

Note that I am no against another solution in general - I don't want to force my idea.

I don't know exactly your proposal. I guess, you want to provide a wrapper type DenseVectorView<Container> that derives from DenseVector, and implements the CRTP interface, i.e., size() and operator[], maybe a few traits classes. This is all fine and would also solve the issues we have here.

But it is not more or less generic than the DenseVectorSpan class, since you make implicit assumptions on the interface of the wrapped Container type, e.g., has a method called size() and operator[] or something. These assumptions are just more explicit in the span implementation, i.e., you have to provide the size and the data accessor/handle to the constructor directly. The deduction guides are just syntactic sugar and could also be discussed to be removed. The span thus allows to use raw arrays, std::vector, std::array, and everything that does not have a std interface, like PETSc containers or armadillo library that uses operator() for element access. It would even make the ScalarVectorView obsolete, since you can also pass a pointer to a scalar and size 1 and it works fine.

I just don't want that we implement yet another span/view/..-like class before we add another one for a tensor generalization (yes, there is also the TensorSpan proposed that generalizes all the views for vectors and matrices that we have currently). The DenseVectorSpan would also be such an intermediate solution, but is has already the same interface. All the additional wrapper classes juts highlight that the implementation in dune-common is somewhat dirty.

An alternative would be to base everything on the iterators. These are quiet generic. Then something like a DenseIteratorRange would solve the problem as well. But then we are also quiet close to the span implementation, just with iterators as data handles and a std::distance for computing the size.

But it is not more or less generic than the DenseVectorSpan class, since you make implicit assumptions on the interface of the wrapped Container type, e.g., has a method called size() and operator[] or something.

In some sense it is more generic, because it does not require to specialize for types as done in the DenseVectorSpan deduction guides. But don't get me wrong: I have the impression that the latter could also be done generically. And if there's a fully generic way to use DenseVectorSpan this would be the preferred solution.

For the time being I'm tempted to add a straight forward implementation of DenseVectorView to Dune::Functions::Impl:: as a temporary hack to fix the issues with activated range checks that make the CI fail for master currently.

as I said, the deduction guides are just syntactic sugar, one maybe should even remove. Instead, one could directly call

auto x = DenseVectorView{inVec.data(), inVec.size()};

and this could even be generalized to

auto x = DenseVectorView{inVec.begin(), inVec.size()};

But this will copy a potentially static size into a runtime variable. Also I cannot see how this can work with proxy containers.

Btw. this is the oldest still open issue in dune-common: core/dune-common#4

mentioned in issue dune-fem/dune-fem#237 (closed)

Fixed by !501 (merged). Notice that we may later want to use DenseTensorSpan and remove the temporary DenseVectorView.

closed