[constraints] Fix contraint computation for hierarchical bases

The algorithm basically interpolates fine DOFs from
coarse basis functions of neighboring elements. This
requires some mechanism to interpolate only hanging fine
DOFs and not all.
Instead of checking if the fine interpolation point is in
within the coarse element, we now check if the interpolated
DOF is associated to a subentity of the intersection.
In most cases this is equivalent. But for hierarchical bases,
the former does probably not lead to the desired result (*).

Remark: (*)
In principle one would expect that the conforming hierarchical
basis consists of the conforming P1-basis plus conforming bubbles.
This is the case for both algorithms (except for some numerical
issues for the old one). However, it turnes out that there are
different possibilities to define the conforming bubbles and
it's not clear which one is desired. Ideally we would like to have:

* The bubbles are lagrangian, i.e. each one is zero at the interpolation
  point of the other ones.
* The bubbles are non-negative.

For hanging nodes you can't have both at once. The old algorithm
guarantees the former, the new one the latter property. I'm undecides
which one is 'correct'.

dune/fufem/constraints/continuityconstraints.hh

@@ -6,8 +6,9 @@
 #include <vector>
 #include <type_traits>
 #include <cmath>
+#include <variant>
 
-#include <dune/geometry/referenceelements.hh>
+#include <dune/functions/functionspacebases/subentitydofs.hh>
 
 #include <dune/fufem/constraints/affineconstraints.hh>
 
@@ -44,26 +45,27 @@ auto makeContinuityConstraints(const Basis& basis)
   // Interpolation weights smaller than the tolerance will be ignored
   const auto tol = 1e-5;
 
-  auto localIndices = [](const auto& tree) {
-    return Dune::range(tree.localIndex(0), tree.localIndex(0)+tree.size());
-  };
-
-  auto localView = basis.localView();
-  auto neighborLocalView = basis.localView();
-
   // Once it was merged in dune-functions the isConstrained vector should be initialized using something like
   // Dune::Functions::istlVectorFactory<char>(std::declval<Basis>().preBasis().containerDescriptor())
   // and the BitVector type should be deduced.
   auto constraints = AffineConstraints<BitVector, MultiIndex, Coefficient>();
-  auto&& interpolation = constraints.interpolation();
   auto&& isConstrained = Dune::Functions::istlVectorBackend(constraints.isConstrained());
-
   isConstrained.resize(basis);
   isConstrained = false;
 
+  auto&& interpolation = constraints.interpolation();
+
   auto localIsConstrained = std::vector<bool>();
   auto interpolationValues = std::vector<Coefficient>();
 
+  auto localIndices = [](const auto& tree) {
+    return Dune::range(tree.localIndex(0), tree.localIndex(0)+tree.size());
+  };
+
+  auto localView = basis.localView();
+  auto neighborLocalView = basis.localView();
+  auto intersectionDOFs = Dune::Functions::subEntityDOFs(basis);
+
   for(const auto& element : elements(basis.gridView()))
   {
     localView.bind(element);
@@ -85,6 +87,8 @@ auto makeContinuityConstraints(const Basis& basis)
       if (neighborElement.level()>=level)
         continue;
 
+      intersectionDOFs.bind(localView, intersection);
+
       neighborLocalView.bind(neighborElement);
 
       Dune::TypeTree::forEachLeafNode(localView.tree(), [&](auto&& node, auto&& treePath) {
@@ -104,15 +108,7 @@ auto makeContinuityConstraints(const Basis& basis)
           neighborBasisValues.resize(k+1);
           auto globalPosition = element.geometry().global(x);
           auto neighborLocalPosition = neighborElement.geometry().local(globalPosition);
-          auto&& neigborRefElement = Dune::referenceElement<double, Element::dimension>(neighborElement.type());
-          if (neigborRefElement.checkInside(neighborLocalPosition))
-            neighborBasis.evaluateFunction(neighborLocalPosition, neighborBasisValues[k]);
-          else
-          {
-            neighborBasisValues[k].resize(neighborBasis.size());
-            for(auto j : Dune::range(neighborBasis.size()))
-              neighborBasisValues[k][j] = 0;
-          }
+          neighborBasis.evaluateFunction(neighborLocalPosition, neighborBasisValues[k]);
           return neighborBasisValues[k++][0];
         };
         node.finiteElement().localInterpolation().interpolate(trackingBasisFunction, interpolationValues);
@@ -121,7 +117,7 @@ auto makeContinuityConstraints(const Basis& basis)
         for(auto j : Dune::range(neighborBasis.size()))
         {
           auto j_global = neighborLocalView.index(neighborNode.localIndex(j));
-          
+
           // neighborBasisFunction_j behaves like the j-th neighbor bases function if
           // interpolate() evaluates it at the same points in the same
           // order as above.
@@ -133,11 +129,16 @@ auto makeContinuityConstraints(const Basis& basis)
 
           for(auto i : Dune::range(node.finiteElement().size()))
           {
+            // Only basis functions associated to this intersection or a subentity
+            // are considered for being interpolated from the neighbor basis.
+            if (not(intersectionDOFs.contains(node.localIndex(i))))
+              continue;
+
             auto i_global = localView.index(node.localIndex(i));
 
             // We need to skip this if dof is already constrained from another element.
             // Otherwise we might get interpolation weights from multiple level neighbors
-            // resulting in doubled interpolation weights after resolving the dependencies
+            // resulting in doubled interpolation weights after resolving the dependencies.
             if (isConstrained[i_global])
               continue;