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

[Simon Praetorius]

Hi Felix,

if you want to access the finiteElement() of a leaf basis node, you have 
to bind it to an element, since, e.g. LagrangeNode uses a 
FiniteElementCache that gives you a different FiniteElement depending on 
the element type. Binding a TypeTree to an element could be done using 
the free function `bindTree(tree, element)` from 
dune/functions/functionspacebases/nodes.hh. But therefore you need a 
mutable tree object. The LocalView gives you only a const& to the 
TypeTree. An option would be to simply create a copy of the tree and 
bind this copy to the element, i.e.

auto tree = lv.tree();
bindTree(tree, coarseElement);
forEachLeafNode(tree, [&](auto const& node, auto&& tp) {
   const auto& fe = node.finiteElement();
   // do something with fe
}/
/

The tree traversal does not require a bound tree, so you could also bind 
the leaf node directly to an element, i.e.

forEachLeafNode(lv.tree(), [&](auto node, auto&& tp) { // node argument 
by value
   node.bind(coarseElement);
   const auto& fe = node.finiteElement();
   // do something with fe
}/
/

This should solve your problem with AlbertaGrid,

Best,
Simon

Am 09.11.18 um 09:07 schrieb 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.
>
>
> _______________________________________________
> dune-functions mailing list
> dune-functions at lists.dune-project.org
> https://lists.dune-project.org/mailman/listinfo/dune-functions
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune-functions/attachments/20181120/da82d4dd/attachment.htm>