Feature/data transfer implementation

src/amdis/DataTransfer.hpp

 #pragma once
 
 #include <memory>
-#include <type_traits>
-#include <utility>
-
-#include <dune/functions/functionspacebases/subspacebasis.hh>
 
 #include <amdis/Output.hpp>
-#include <amdis/utility/TreeData.hpp>
-#include <amdis/utility/Visitor.hpp>
 
 namespace AMDiS
 {
@@ -17,17 +11,6 @@ namespace AMDiS
     INTERPOLATE  = 1
   } DataTransferOperation;
 
-  template <class Node, class RangeType>
-  class NodeDataTransfer
-  {
-  public:
-    template <class Basis, class Container>
-    void preAdapt(Basis const& basis, Container const& coeff, bool mightCoarsen) {}
-
-    template <class Basis, class Container>
-    void postAdapt(Basis const& basis, Container& coeff, bool refined) const {}
-  };
-
 
   template <class Container>
   class DataTransferInterface
@@ -40,7 +23,7 @@ namespace AMDiS
     virtual void preAdapt(Container const& container, bool mightCoarsen) = 0;
 
     /// Interpolate data to new grid after grid adaption
-    virtual void postAdapt(Container& container, bool refined) const = 0;
+    virtual void postAdapt(Container& container, bool refined) = 0;
   };
 
 
@@ -51,9 +34,9 @@ namespace AMDiS
       : public DataTransferInterface<Container>
   {
   public:
-    virtual void preAdapt(Container const& container, bool) override {}
+    void preAdapt(Container const& container, bool) override {}
 
-    virtual void postAdapt(Container& container, bool) const override
+    void postAdapt(Container& container, bool) override
     {
       container.compress();
     }
@@ -61,45 +44,7 @@ namespace AMDiS
 
 
   template <class Container, class Basis>
-  class DataTransfer
-      : public DataTransferInterface<Container>
-  {
-    template <class Node>
-    using NDT = NodeDataTransfer<std::decay_t<Node>, typename Container::value_type>;
-
-  public:
-    DataTransfer(Basis const& basis)
-      : basis_(&basis)
-    {}
-
-    /// Calls \ref NodeDataTransfer::preAdapt() on each basis node
-    virtual void preAdapt(Container const& container, bool mightCoarsen) override
-    {
-      nodeDataTransfer_.init(*basis_);
-      AMDiS::forEachLeafNode_(basis_->localView().tree(), [&](auto const& node, auto const& treePath)
-      {
-        auto subBasis = Dune::Functions::subspaceBasis(*basis_, treePath);
-        nodeDataTransfer_[node].preAdapt(subBasis, container.vector(), mightCoarsen);
-      });
-    }
-
-    /// Calls \ref NodeDataTransfer::postAdapt() on each basis node after compressing the
-    /// Container the dimension of the basis
-    virtual void postAdapt(Container& container, bool refined) const override
-    {
-      container.compress();
-      AMDiS::forEachLeafNode_(basis_->localView().tree(), [&](auto const& node, auto const& treePath)
-      {
-        auto subBasis = Dune::Functions::subspaceBasis(*basis_, treePath);
-        nodeDataTransfer_[node].postAdapt(subBasis, container.vector(), refined);
-      });
-    }
-
-  private:
-    Basis const* basis_;
-    TreeData<Basis, NDT, true> nodeDataTransfer_;
-  };
-
+  class DataTransfer;
 
   /// Factory to create DataTransfer objects based on the \ref DataTransferOperation
   template <class Container>
@@ -116,7 +61,7 @@ namespace AMDiS
         case NO_OPERATION:
           return std::make_unique<NoDataTransfer<Container>>();
         case INTERPOLATE:
-          return std::make_unique<DataTransfer<Container, Basis>>(basis);
+          return std::make_unique< DataTransfer<Container, Basis> >(basis);
         default:
           error_exit("Invalid data transfer\n");
           return nullptr; // avoid warnings
@@ -125,3 +70,5 @@ namespace AMDiS
   };
 
 } // end namespace AMDiS
+
+#include "DataTransfer.inc.hpp"

src/amdis/DataTransfer.inc.hpp

+#pragma once
+
+#include <cmath>
+#include <functional>
+#include <limits>
+#include <map>
+#include <memory>
+#include <numeric>
+#include <type_traits>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/hash.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/common/rangegenerators.hh>
+
+#include <dune/functions/common/functionfromcallable.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <amdis/Output.hpp>
+#include <amdis/utility/ConcurrentCache.hpp>
+#include <amdis/utility/TreeContainer.hpp>
+#include <amdis/utility/Visitor.hpp>
+
+namespace AMDiS
+{
+  template <class Node, class Container, class Basis>
+  class NodeDataTransfer {};
+
+
+  /** Data Transfer implementation for a single grid
+   *  Handles computations related to the geometric information of the grid and passes that to the
+   *    underlying NodeDataTransfer classes
+   */
+  template <class Container, class Basis>
+  class DataTransfer
+      : public DataTransferInterface<Container>
+  {
+    using LocalView = typename Basis::LocalView;
+    using Tree = typename LocalView::Tree;
+    using GridView = typename Basis::GridView;
+    using Grid = typename GridView::Grid;
+    using Mapper = Dune::LeafMultipleCodimMultipleGeomTypeMapper<Grid>;
+
+    using IdSet = typename Grid::LocalIdSet;
+    using IdType = typename IdSet::IdType;
+
+    using Element = typename GridView::template Codim<0>::Entity;
+    using Geometry = typename Element::Geometry;
+    using LocalCoordinate = typename Geometry::LocalCoordinate;
+    using BoolCoordPair = std::pair<bool, LocalCoordinate>;
+
+    // Hash function for cache container
+    struct Hasher
+    {
+      std::size_t operator()(LocalCoordinate const& coord) const
+      {
+        std::size_t seed = 0;
+        for (std::size_t i = 0; i < coord.size(); ++i)
+          Dune::hash_combine(seed, coord[i]);
+        return seed;
+      }
+    };
+
+    using CacheImp = std::unordered_map<LocalCoordinate, BoolCoordPair, Hasher>;
+    using ChildCache = ConcurrentCache<LocalCoordinate, BoolCoordPair, ConsecutivePolicy, CacheImp>;
+
+    constexpr static int d = Geometry::coorddimension;
+    using ctype = typename Geometry::ctype;
+
+    // Returns the Node's NodeDataTransfer
+    struct NDT
+    {
+      template<class Node>
+      auto operator()(const Node& node)
+      {
+        using NDT = NodeDataTransfer<Node, Container, Basis>;
+        return NDT{};
+      }
+    };
+
+    using NodeDataTransferContainer = std::decay_t<decltype(
+      makeTreeContainer<Tree, NDT>(std::declval<const Tree&>(), NDT()))>;
+
+    // Returns the Node's NodeElementData
+    struct NodeElementData
+    {
+      template<class Node>
+      auto operator()(const Node& node)
+      {
+        using NED = typename NodeDataTransfer<Node, Container, Basis>
+          ::NodeElementData;
+        return NED{};
+      }
+    };
+
+    using ElementData = std::decay_t<decltype(
+      makeTreeContainer<Tree, NodeElementData>(std::declval<const Tree&>(), NodeElementData()))>;
+
+  public:
+    // Container with data that persists during grid adaptation
+    using PersistentContainer = std::map<IdType, ElementData>;
+
+  public:
+    DataTransfer(Basis const& basis)
+      : basis_(&basis)
+      , mapper_(basis.gridView().grid(), Dune::mcmgElementLayout())
+      , finished_(mapper_.size(), false)
+      , nodeDataTransfer_(makeTreeContainer<Tree, NDT>(basis_->localView().tree(), NDT()))
+    {}
+
+    /** Saves data contained in coeff in the PersistentContainer
+     *  To be called after grid.preAdapt() and before grid.adapt()
+     */
+    void preAdapt(Container const& coeff, bool mightCoarsen) override;
+
+    /** Unpacks data from the PersistentContainer
+     *  To be called after grid.adapt() and before grid.postAdapt()
+     */
+    void postAdapt(Container& coeff, bool refined) override;
+
+  private:
+    Basis const* basis_;
+    PersistentContainer persistentContainer_;
+    Mapper mapper_;
+    std::vector<bool> finished_;
+    NodeDataTransferContainer nodeDataTransfer_;
+  };
+
+  template <class Container, class Basis>
+  void DataTransfer<Container, Basis>::preAdapt(Container const& coeff, bool mightCoarsen)
+  {
+    auto gv = basis_->gridView();
+    auto lv = basis_->localView();
+    auto const& idSet = gv.grid().localIdSet();
+
+    forEachLeafNode_(lv.tree(), [&](auto const& node, auto const& tp) {
+      nodeDataTransfer_[tp].preAdaptInit(lv, coeff, node);
+    });
+
+    // Make persistent DoF container
+    persistentContainer_.clear();
+    for (const auto& e : elements(gv))
+    {
+      auto it = persistentContainer_.emplace(idSet.id(e),
+        makeTreeContainer<Tree, NodeElementData>(lv.tree(), NodeElementData()));
+
+      lv.bind(e);
+      auto& treeContainer = it.first->second;
+      forEachLeafNode_(lv.tree(), [&](auto const& node, auto const& tp) {
+        nodeDataTransfer_[tp].cacheLocal(treeContainer[tp]);
+      });
+    }
+
+    // Interpolate from possibly vanishing elements
+    if (mightCoarsen) {
+      auto maxLevel = gv.grid().maxLevel();
+      ctype const checkInsideTolerance = std::sqrt(std::numeric_limits<ctype>::epsilon());
+      for (auto const& e : elements(gv))
+      {
+        auto father = e;
+        while (father.mightVanish() && father.hasFather())
+        {
+          father = father.father();
+          auto it = persistentContainer_.emplace(idSet.id(father),
+            makeTreeContainer<Tree, NodeElementData>(lv.tree(), NodeElementData()));
+          if (it.second) {
+            auto& treeContainer = it.first->second;
+            auto geo = father.geometry();
+            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) {
+              if (hIt->isLeaf()) {
+                auto const& child = *hIt;
+                auto search = persistentContainer_.find(idSet.id(child));
+                assert(search != persistentContainer_.end());
+                auto const& childContainer = search->second;
+                lv.bind(child);
+
+                auto const& childGeo = child.geometry();
+                auto const& ref = Dune::referenceElement(childGeo);
+                auto const& refTypeId = ref.type().id();
+
+                // 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::template checkInside<ctype,d>
+                                    (refTypeId, d, local, checkInsideTolerance);
+                  return BoolCoordPair(isInside, std::move(local));
+                };
+
+                // TODO(FM): Disable for single-node basis
+                ChildCache childCache;
+                auto xInChildCached = [&](LocalCoordinate const& x) -> BoolCoordPair {
+                  return childCache.get(x, [&](LocalCoordinate const& x) { return xInChild(x); });
+                };
+
+                restrictLocalCompleted = true;
+                forEachLeafNode_(lv.tree(), [&](auto const& node, auto const& tp) {
+                  restrictLocalCompleted &=
+                    nodeDataTransfer_[tp].restrictLocal(father, treeContainer[tp], xInChildCached,
+                                                        childContainer[tp], init);
+                });
+                init = false;
+              }
+            }
+            // test if restrictLocal was completed on all nodes
+            assert(restrictLocalCompleted);
+          }
+        }
+      }
+    }
+  }
+
+  template <class Container, class Basis>
+  void DataTransfer<Container, Basis>::postAdapt(Container& coeff, bool refined)
+  {
+    coeff.resize(*basis_);
+    auto gv = basis_->gridView();
+    auto lv = basis_->localView();
+    auto const& idSet = gv.grid().localIdSet();
+    forEachLeafNode_(lv.tree(), [&](auto const& node, auto const& tp) {
+      nodeDataTransfer_[tp].postAdaptInit(lv, coeff, node);
+    });
+
+    mapper_.update();
+    finished_.clear();
+    finished_.resize(mapper_.size(), false);
+    for (const auto& e : elements(gv))
+    {
+      auto index = mapper_.index(e);
+      if (finished_[index])
+        continue;
+
+      auto it = persistentContainer_.find(idSet.id(e));
+
+      // Data already exists and no interpolation is required
+      if (it != persistentContainer_.end()) {
+        lv.bind(e);
+        auto const& treeContainer = it->second;
+        forEachLeafNode_(lv.tree(), [&](auto const& node, auto const& tp) {
+          nodeDataTransfer_[tp].copyLocal(treeContainer[tp]);
+        });
+        finished_[index] = true;
+      }
+      // Data needs to be interpolated
+      else {
+        auto father = e;
+        while (father.hasFather() && father.isNew())
+          father = father.father();
+
+        auto maxLevel = gv.grid().maxLevel();
+        auto fatherGeo = father.geometry();
+        bool init = true; // init flag for first call on new father element
+
+        auto it = persistentContainer_.find(idSet.id(father));
+        assert(it != persistentContainer_.end());
+        auto const& treeContainer = it->second;
+
+        auto hItEnd = father.hend(maxLevel);
+        for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt) {
+          if (hIt->isLeaf()) {
+            auto const& child = *hIt;
+            lv.bind(child);
+
+            // coordinate transform from child to father element
+            auto xInFather = [&fatherGeo, childGeo = child.geometry()]
+              (LocalCoordinate const& x) -> LocalCoordinate
+            {
+              return fatherGeo.local(childGeo.global(x));
+            };
+
+            forEachLeafNode_(lv.tree(), [&](auto const& node, auto const& tp) {
+              nodeDataTransfer_[tp].prolongLocal(father, treeContainer[tp], xInFather, init);
+            });
+
+            finished_[mapper_.index(child)] = true;
+            init = false;
+          }
+        }
+      }
+    }
+  }
+
+  /** Element-local data transfer on a single leaf node of the basis tree
+   *  Handles computations related to the finite element basis node
+   */
+  template <class GV, int k, class TP, class Container, class Basis>
+  class NodeDataTransfer<Dune::Functions::LagrangeNode<GV,k,TP>, Container, Basis>
+  {
+    using T = typename Container::value_type;
+    using LocalView = typename Basis::LocalView;
+    using Element = typename GV::template Codim<0>::Entity;
+
+    using Node = typename Dune::Functions::LagrangeNode<GV,k,TP>;
+    using LocalBasis = typename Node::FiniteElement::Traits::LocalBasisType;
+    using LBRangeType = typename LocalBasis::Traits::RangeType;
+    using LocalInterpolation = typename Node::FiniteElement::Traits::LocalBasisType;
+    using LIDomainType = typename LocalInterpolation::Traits::DomainType;
+    using LIRangeType = typename LocalInterpolation::Traits::RangeType;
+
+  public:
+    using NodeElementData = std::vector<T>;
+
+  public:
+    NodeDataTransfer() = default;
+
+    /// To be called once before cacheLocal/restrictLocal are called within the preAdapt step
+    void preAdaptInit(LocalView const& lv, Container const& coeff, Node const& node)
+    {
+      lv_ = &lv;
+      node_ = &node;
+      auto nodeCopy = node;
+      fatherNode_ = std::make_unique<Node>(std::move(nodeCopy));
+      constCoeff_ = &coeff;
+    }
+
+    /// To be called once before copyLocal/prolongLocal are called within the postAdapt step
+    void postAdaptInit(LocalView const& lv, Container& coeff, Node const& node)
+    {
+      lv_ = &lv;
+      node_ = &node;
+      auto nodeCopy = node;
+      fatherNode_ = std::make_unique<Node>(std::move(nodeCopy));
+      coeff_ = &coeff;
+    }
+
+    /// Cache data on the element bound to node_
+    void cacheLocal(NodeElementData& dofs) const
+    {
+      auto const& fe = node_->finiteElement();
+      auto feSize = fe.size();
+      dofs.resize(feSize);
+      for (std::size_t i = 0; i < feSize; ++i)
+        dofs[i] = (*constCoeff_)[lv_->index(node_->localIndex(i))];
+    }
+
+    /** Evaluate data on the child element bound to node_ and interpolate onto father entity
+     *    using the coordinate transformation trafo from father to child
+     */
+    template <class Trafo>
+    bool restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& trafo,
+                       NodeElementData const& childDOFs, bool init);
+
+    /// Copy already existing data to element bound to node_
+    void copyLocal(NodeElementData const& dofs) const
+    {
+      for (std::size_t i = 0; i < dofs.size(); ++i)
+        (*coeff_)[lv_->index(node_->localIndex(i))] = dofs[i];
+    }
+
+    /** Interpolate data from father onto the child element bound to node_ using the transformation
+     *    trafo from child to father
+     */
+    template <class Trafo>
+    void prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
+                      Trafo const& trafo, bool init);
+
+  private:
+    LocalView const* lv_;
+    Node const* node_;
+    std::unique_ptr<Node> fatherNode_;
+    Container const* constCoeff_;
+    Container* coeff_;
+    std::vector<bool> finishedDOFs_;
+    NodeElementData fatherDOFsTemp_;
+  };
+
+  template <class GV, int k, class TP, class Container, class Basis>
+    template <class Trafo>
+  bool NodeDataTransfer<Dune::Functions::LagrangeNode<GV,k,TP>, Container, Basis>::
+    restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& trafo,
+                  NodeElementData const& childDOFs, bool init)
+  {
+    auto& fatherNode = *fatherNode_;
+    std::size_t currentDOF = 0;
+    if (init)
+    {
+      // TODO(FM): This is UB, replace with FE cache for father
+      bindTree(fatherNode, father);
+    }
+    auto const& childNode = *node_;
+    auto const& childFE = childNode.finiteElement();
+    auto const& fatherFE = fatherNode.finiteElement();
+    if (init)
+    {
+      finishedDOFs_.assign(fatherFE.size(), false);
+      fatherDOFsTemp_.assign(fatherFE.size(), 0);
+    }
+
+    auto evalLeaf = [&](LIDomainType const& x) -> LIRangeType {
+      if (!finishedDOFs_[currentDOF])
+      {
+        auto const& insideLocal = trafo(x);
+        bool isInside = insideLocal.first;
+        if (isInside)
+        {
+          auto const& local = insideLocal.second;
+          thread_local std::vector<LBRangeType> shapeValues;
+          childFE.localBasis().evaluateFunction(local, shapeValues);
+
+          assert(childDOFs.size() == shapeValues.size());
+
+          LIRangeType y(0);
+          for (std::size_t i = 0; i < shapeValues.size(); ++i)
+            y += shapeValues[i] * childDOFs[i];
+
+          fatherDOFsTemp_[currentDOF] = y;
+          finishedDOFs_[currentDOF++] = true;
+          return y;
+        }
+      }
+      return fatherDOFsTemp_[currentDOF++];
+    };
+
+    using FFC
+      = Dune::Functions::FunctionFromCallable<LIRangeType(LIDomainType), decltype(evalLeaf)>;
+    fatherFE.localInterpolation().interpolate(FFC(evalLeaf), fatherDOFs);
+
+    // Return true if all father DOFs have been evaluated
+    return std::accumulate(finishedDOFs_.begin(), finishedDOFs_.end(), true,
+                           std::logical_and<bool>());
+  }
+
+  template <class GV, int k, class TP, class Container, class Basis>
+    template <class Trafo>
+  void NodeDataTransfer<Dune::Functions::LagrangeNode<GV,k,TP>, Container, Basis>::
+    prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
+                 Trafo const& trafo, bool init)
+  {
+    auto& fatherNode = *fatherNode_;
+    if (init)
+    {
+      // TODO(FM): This is UB, replace with FE cache for father
+      bindTree(fatherNode, father);
+    }
+    auto const& childNode = *node_;
+
+    // evaluate father in child coordinate x
+    auto evalFather = [&](LIDomainType const& x) -> LIRangeType
+    {
+      thread_local std::vector<LBRangeType> shapeValues;
+      fatherNode.finiteElement().localBasis().evaluateFunction(trafo(x), shapeValues);
+      assert(shapeValues.size() == fatherDOFs.size());
+
+      LIRangeType y(0);
+      for (std::size_t i = 0; i < shapeValues.size(); ++i)
+        y += shapeValues[i] * fatherDOFs[i];
+
+      return y;
+    };
+
+    auto const& childFE = childNode.finiteElement();
+    thread_local std::vector<T> childDOFs;
+    using FFC = Dune::Functions
+      ::FunctionFromCallable<LIRangeType(LIDomainType), decltype(evalFather)>;
+    childFE.localInterpolation().interpolate(FFC(evalFather), childDOFs);
+
+    for (std::size_t i = 0; i < childDOFs.size(); ++i)
+      (*coeff_)[lv_->index(childNode.localIndex(i))] = childDOFs[i];
+  }
+
+} // end namespace AMDiS

src/amdis/common/Apply.hpp

+#pragma once
+
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+#include <amdis/common/Mpl.hpp>
+
+namespace AMDiS
+{
+  namespace Tools
+  {
+    namespace Impl_
+    {
+      template <class F, class Tuple, std::size_t... I>
+      constexpr decltype(auto) apply_impl(F&& f, Tuple&& t, std::index_sequence<I...>)
+      {
+        return f(std::get<I>(std::forward<Tuple>(t))...);
+      }
+
+      template <class F, std::size_t I0, std::size_t... I>