Merge branch 'fix/datatransfer_geoInFather' into 'master'

Fix issues with InterpolationDataTransfer on curved grids

See merge request amdis/amdis!90

amdis/datatransfers/InterpolationDataTransfer.inc.hpp

@@ -50,7 +50,8 @@ preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
 {
   GridView gv = basis.gridView();
   LocalView lv = basis.localView();
-  auto const& idSet = gv.grid().localIdSet();
+  auto const& grid = gv.grid();
+  auto const& idSet = grid.localIdSet();
 
   Traversal::forEachLeafNode(lv.tree(), [&](auto const& node, auto const& tp) {
     nodeDataTransfer_[tp].preAdaptInit(lv, coeff, node);
@@ -73,9 +74,12 @@ preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
     return;
 
   // Interpolate from possibly vanishing elements
-  auto maxLevel = gv.grid().maxLevel();
+  auto maxLevel = grid.maxLevel();
   using std::sqrt;
   typename Grid::ctype const checkInsideTolerance = sqrt(std::numeric_limits<typename Grid::ctype>::epsilon());
+  using Seed = typename Grid::template Codim<0>::EntitySeed;
+  std::vector<Seed> seeds(maxLevel+1);
+
   for (auto const& e : entitySet(basis))
   {
     auto father = e;
@@ -91,7 +95,10 @@ preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
       bool init = true; // init flag for first call on new father element
       bool restrictLocalCompleted = false;
       auto hItEnd = father.hend(maxLevel);
-      for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt) {
+      for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt)
+      {
+        seeds[hIt->level()] = hIt->seed();  // Save element in hierarchy to access geometryInFather at each step later
+
         if (!hIt->isLeaf())
           continue;
 
@@ -111,12 +118,21 @@ preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
         // Transfers input father-local point x into child-local point y
         // Returns false if x is not inside the child
         auto xInChild = [&](LocalCoordinate const& x) -> BoolCoordPair {
-          LocalCoordinate local = childGeo.local(geo.global(x));
-          // TODO(FM): Using an implementation detail as workaround for insufficient
-          //   tolerance, see https://gitlab.dune-project.org/core/dune-grid/issues/84
-          bool isInside = Dune::Geo::Impl::checkInside(refTypeId, Geometry::mydimension, local, checkInsideTolerance);
-          return BoolCoordPair(isInside, std::move(local));
+          LocalCoordinate y = x;
+          for (std::size_t i = father.level()+1; i <= child.level(); ++i)
+          {
+            auto currentElement = grid.entity(seeds[i]);
+            auto const& geoInFather = currentElement.geometryInFather();
+            y = geoInFather.local(y);
+            // TODO(FM): Using an implementation detail as workaround for insufficient
+            //   tolerance, see https://gitlab.dune-project.org/core/dune-grid/issues/84
+            bool isInside = Dune::Geo::Impl::checkInside(Dune::referenceElement(geoInFather).type().id(), Geometry::mydimension, y, checkInsideTolerance);
+            if (!isInside)
+              return BoolCoordPair(false, std::move(y));
+          }
+          return BoolCoordPair(true, std::move(y));
         };
+        // Cache result of xInChild for subsequent calls from other basis nodes
         // TODO(FM): Disable for single-node basis
         ChildCache childCache;
         auto xInChildCached = [&](LocalCoordinate const& x) -> BoolCoordPair {
@@ -201,10 +217,16 @@ void InterpolationDataTransfer<B,C>::adapt(B const& basis, C& coeff)
       lv.bind(child);
 
       // coordinate transform from child to father element
-      auto xInFather = [&fatherGeo, childGeo = child.geometry()]
-        (LocalCoordinate const& x) -> LocalCoordinate
+      auto xInFather = [&](LocalCoordinate const& x) -> LocalCoordinate
       {
-        return fatherGeo.local(childGeo.global(x));
+        auto y = x;
+        auto currentElement = child;
+        while (currentElement.level() != father.level())
+        {
+          y = currentElement.geometryInFather().global(y);
+          currentElement = currentElement.father();
+        }
+        return y;
       };
 
       Traversal::forEachLeafNode(lv.tree(), [&](auto const& /*node*/, auto const& tp) {