[dune-functions] Binding to non-leaf entities in dune-functions 2.6

[Felix Mueller]

Dear dune-functions developers,

I am working on some code for data transfer during grid adaptation. In this
I have a data vector that contains the degrees of freedom (DOF) of the
grid's finite elements, modeled by a vector indexed with the index set of a
dune-functions basis.
Since I am losing geometric grid information of elements that disappear
during a call to adapt(), I want to compute the DOFs of coarse grid
entities before calling adapt(). To do this I require the finite element
type of those entities, as well as the indices corresponding to the child
entites of the coarse entity.

The following code illustrates the process:

/auto lv = basis.localView();
// loop over all leaf elements
for (auto leafElement : leafElements) {
  lv.bind(leafElement);
  const auto& node = lv.tree();
  const auto& fe = node.finiteElement();
  auto feSize = fe.size();

  for (std::size_t i = 0; i < feSize; ++i) {
    auto index = lv.index(node.localIndex(i);
    // save vector[index] for later use
    // ...
  }
}

lv.bind(coarseElement);
const auto& node = lv.tree();
const auto& fe = node.finiteElement();
// do something with fe
// .../

Note that /coarseElement /does not access the basis's index set here.

This process works for /Dune::UGGrid<2> /and
/Dune::Functions::LagrangeBasis<GridView, k> /with /k/ of any order and
Dune::AlbertaGrid<2,2> with /Dune::Functions::LagrangeBasis<GridView, 1>/.

However it fails for Dune::AlbertaGrid<2,2> with
/Dune::Functions::LagrangeBasis<GridView, k> /with/ k ≠ 1/ at:
DUNE/GRID/ALBERTAGRID/INDEXSETS.HH:482:/ Assertion `(subIndex >= 0) &&
(subIndex < size( codim ))' failed./
Deactivating the assertion prevents the crash and makes my code run just
fine.

The issue above seems to be caused by
DUNE/FUNCTIONS/FUNCTIONSPACEBASES/DEFAULTLOCALVIEW.HH:72:
void bind(const Element& e)
{
element_ = e;
bindTree(tree_, element_);
nodeIndexSet_.bind(tree_);
indices_.resize(size());
Hybrid::ifElse(
Std::is_detected<hasIndices,NodeIndexSet>{},
[&](auto id) {
id(nodeIndexSet_).indices(indices_.begin());
},
[&](auto id) {
for (size_type i = 0 ; i < this->size() ; ++i)
indices_[i] = id(nodeIndexSet_).index(i);
});
}
This tries to access the index set on /coarseElement/, which is not defined
and not needed within my code either.
DUNE-FUNCTIONS 2.5 had a seperate binding mechanism for the local view and
index sets, providing only the method
/void bind(const Element& e)
{
  element_ = e;
  bindTree(tree_, element_);
}/
which is completely sufficient for my purpose.

_My question now is:_
- Is this behaviour intended or is this a bug in either DEFAULTLOCALVIEW.HH
or ALBERTAGRID?
- If it is intended, is there a way to only bind the local view without
binding the index set, like in v2.5?

Thanks for your attention and helpful answers!

Felix Müller

PS: You will find a small file reproducing the error attached. A filename
for the AlbertaReader must be provided as argument to the main method. I
have used DUNE-GRID/DOC/GRIDS/AMC/GRID-2-2.AMC as grid file for my tests.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune-functions/attachments/20181109/0fdbb820/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: AlbertaIndexSetBug.cc
Type: text/x-c
Size: 1655 bytes
Desc: not available
URL: <https://lists.dune-project.org/pipermail/dune-functions/attachments/20181109/0fdbb820/attachment.bin>