InterpolationDataTransfer.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/common/version.hh>

#include <dune/grid/common/geometry.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/rangegenerators.hh>

#include <amdis/Output.hpp>
#include <amdis/common/ConcurrentCache.hpp>
#include <amdis/functions/EntitySet.hpp>
#include <amdis/operations/Assigner.hpp>
#include <amdis/typetree/Traversal.hpp>

namespace AMDiS {
namespace Impl {

// Hash function for cache container
struct CoordHasher
{
  template <class LocalCoord>
  std::size_t operator()(LocalCoord 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;
  }
};

} // end namespace Impl


template <class B, class C>
void InterpolationDataTransfer<B,C>::
preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
{
  GridView gv = basis.gridView();
  LocalView lv = basis.localView();
  auto const& grid = gv.grid();
  auto const& idSet = grid.localIdSet();

  nodeDataTransfer_.resize(lv.tree());
  Traversal::forEachLeafNode(lv.tree(), [&](auto const& node, auto const& tp) {
    nodeDataTransfer_[tp].preAdaptInit(lv, coeff, node);
  });

  // Make persistent DoF container
  persistentContainer_.clear(); // Redundant if postAdapt was correctly called last cycle
  for (const auto& e : elements(gv))
  {
    auto it = persistentContainer_.emplace(idSet.id(e),
      Dune::TypeTree::makeTreeContainer(lv.tree(), [](auto&& node) { return NodeElementData<TYPEOF(node)>{}; }));

    lv.bind(e);
    auto& treeContainer = it.first->second;
    Traversal::forEachLeafNode(lv.tree(), [&](auto const& /*node*/, auto const& tp) {
      nodeDataTransfer_[tp].cacheLocal(treeContainer[tp]);
    });
  }

  if (!mightCoarsen)
    return;

  // Interpolate from possibly vanishing elements
  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;
    while (father.mightVanish() && father.hasFather())
    {
      father = father.father();
      auto it = persistentContainer_.emplace(idSet.id(father),
        Dune::TypeTree::makeTreeContainer(lv.tree(), [](auto&& node) { return NodeElementData<TYPEOF(node)>{}; }));
      if (!it.second)
        continue;

      auto& treeContainer = it.first->second;
      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)
      {
        seeds[hIt->level()] = hIt->seed();  // Save element in hierarchy to access geometryInFather at each step later

        if (!hIt->isLeaf())
          continue;

        auto const& child = *hIt;
        auto search = persistentContainer_.find(idSet.id(child));
        assert(search != persistentContainer_.end());
        auto const& childContainer = search->second;
        lv.bind(child);

        using BoolCoordPair = std::pair<bool, LocalCoordinate>;
        using CacheImp = std::unordered_map<LocalCoordinate, BoolCoordPair, Impl::CoordHasher>;
        using ChildCache = ConcurrentCache<LocalCoordinate, BoolCoordPair, ConsecutivePolicy, CacheImp>;

        // 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 y = x;
          for (int 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 {
          return childCache.get(x, [&](LocalCoordinate const& x) { return xInChild(x); });
        };

        restrictLocalCompleted = true;
        Traversal::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);

    } // end while (father.mightVanish)
  } // end for (elements)
}


template <class B, class C>
void InterpolationDataTransfer<B,C>::adapt(B const& basis, C& coeff)
{
  coeff.resize(basis);

  // No data was saved before adapting the grid, make
  // sure to call DataTransfer::preAdapt before calling adapt() on the grid
  if (persistentContainer_.empty())
    return;

  GridView gv = basis.gridView();
  LocalView lv = basis.localView();
  auto const& idSet = gv.grid().localIdSet();
  Traversal::forEachLeafNode(lv.tree(), [&](auto const& node, auto const& tp) {
    nodeDataTransfer_[tp].adaptInit(lv, coeff, node);
  });

  using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
  Mapper mapper{gv, Dune::mcmgElementLayout()};

  std::vector<bool> finished(mapper.size(), false);
  for (const auto& e : entitySet(basis))
  {
    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;
      Traversal::forEachLeafNode(lv.tree(), [&](auto const& /*node*/, auto const& tp) {
        nodeDataTransfer_[tp].copyLocal(treeContainer[tp]);
      });
      finished[index] = true;
      continue;
    }

    // Data needs to be interpolated
    auto father = e;
    while (father.hasFather() && father.isNew())
      father = father.father();

    auto maxLevel = gv.grid().maxLevel();
    bool init = true; // init flag for first call on new father element

    auto father_it = persistentContainer_.find(idSet.id(father));
    assert(father_it != persistentContainer_.end());
    auto const& treeContainer = father_it->second;

    auto hItEnd = father.hend(maxLevel);
    for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt) {
      if (!hIt->isLeaf())
        continue;

      auto const& child = *hIt;
      lv.bind(child);

      // coordinate transform from child to father element
      auto xInFather = [&](LocalCoordinate const& x) -> LocalCoordinate
      {
        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) {
        nodeDataTransfer_[tp].prolongLocal(father, treeContainer[tp], xInFather, init);
      });

      finished[mapper.index(child)] = true;
      init = false;
    }
  } // end for (elements)

  coeff.finish();
}


template <class B, class C>
void  InterpolationDataTransfer<B,C>::postAdapt(C&)
{
  persistentContainer_.clear();
}


template <class Node, class Container, class Basis>
class NodeDataTransfer
{
  using T = typename Container::value_type;
  using LocalView = typename Basis::LocalView;
  using Element = typename Node::Element;

  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;

  void preAdaptInit(LocalView const& lv, Container const& coeff, Node const& node)
  {
    lv_ = &lv;
    node_ = &node;
    fatherNode_ = std::make_shared<Node>(node);
    constCoeff_ = &coeff;
  }


  // [[expects: preAdaptInit to be called before]]
  void cacheLocal(NodeElementData& dofs) const
  {
    constCoeff_->gather(*lv_, *node_, dofs);
  }

  // [[expects: preAdaptInit to be called before]]
  template <class Trafo>
  bool restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& trafo,
                      NodeElementData const& childDOFs, bool init);


  void adaptInit(LocalView const& lv, Container& coeff, Node const& node)
  {
    lv_ = &lv;
    node_ = &node;
    fatherNode_ = std::make_shared<Node>(node);
    coeff_ = &coeff;
  }

  // [[expects: adaptInit to be called before]]
  void copyLocal(NodeElementData const& dofs) const
  {
    coeff_->scatter(*lv_, *node_, dofs, Assigner::assign{});
  }

  // [[expects: adaptInit to be called before]]
  template <class Trafo>
  void prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
                    Trafo const& trafo, bool init);

private:
  LocalView const* lv_ = nullptr;
  Node const* node_ = nullptr;
  std::shared_ptr<Node> fatherNode_;
  Container const* constCoeff_ = nullptr;
  Container* coeff_ = nullptr;
  std::vector<bool> finishedDOFs_;
  NodeElementData fatherDOFsTemp_;
};


template <class N, class C, class B>
  template <class Trafo>
bool NodeDataTransfer<N,C,B>::
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, see https://gitlab.com/amdis/amdis/-/issues/16 (case 2)
    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] = T(y);
        finishedDOFs_[currentDOF++] = true;
        return y;
      }
    }
    return fatherDOFsTemp_[currentDOF++];
  };

  fatherFE.localInterpolation().interpolate(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 N, class C, class B>
  template <class Trafo>
void NodeDataTransfer<N,C,B>::
prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
             Trafo const& trafo, bool init)
{
  auto& fatherNode = *fatherNode_;
  if (init)
  {
    // TODO(FM): This is UB, see https://gitlab.com/amdis/amdis/-/issues/16 (case 1)
    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;
  childFE.localInterpolation().interpolate(evalFather, childDOFs);

  coeff_->scatter(*lv_, childNode, childDOFs, Assigner::assign{});
}

} // end namespace AMDiS