SimpleInterpolator.hpp

#pragma once
 
#include <vector>
 
#include <amdis/common/FakeContainer.hpp>
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/functions/EntitySet.hpp>
#include <amdis/functions/HierarchicNodeToRangeMap.hpp>
#include <amdis/functions/NodeIndices.hpp>
#include <amdis/linearalgebra/Traits.hpp>
#include <amdis/operations/Assigner.hpp>
#include <amdis/typetree/Traversal.hpp>
 
#include <dune/common/typetree/childaccess.hh>
#include <dune/common/typetree/treepath.hh>
 
namespace AMDiS
{
  namespace tag
  {
    struct assign {};
  }
 
  template <class Basis,
            class TreePath = Dune::TypeTree::TreePath<>>
  struct SimpleInterpolator
  {
    template <class Coeff, class GridFct, class BitVector>
    void operator()(Coeff& coeff, GridFct const& gf, BitVector const& bitVec) const
    {
      // Obtain a local view of the gridFunction
      auto lf = localFunction(gf);
      auto localView = basis_.localView();
 
      std::vector<typename Coeff::value_type> localCoeff;
 
      for (const auto& e : entitySet(basis_))
      {
        localView.bind(e);
        lf.bind(e);
 
        auto&& subTree = Dune::TypeTree::child(localView.tree(),treePath_);
        Traversal::forEachLeafNode(subTree, [&](auto const& node, auto const& tp)
        {
          auto bitVecRange = mappedRangeView(Dune::range(node.size()), [&](std::size_t i) -> bool {
            return bitVec[localView.index(node.localIndex(i))];
          });
 
          // check whether there is a local contribution
          std::vector<bool> mask(bitVecRange.begin(), bitVecRange.end());
          if (std::all_of(mask.begin(), mask.end(), [](bool m) { return !m; }))
            return;
 
          // compute local interpolation coefficients
          node.finiteElement().localInterpolation().interpolate(HierarchicNodeWrapper{tp,lf}, localCoeff);
 
          // assign local coefficients to global vector
          coeff.scatter(localView, node, localCoeff, mask, Assigner::assign{});
        });
 
        lf.unbind();
      }
      coeff.finish();
    }
 
    template <class Coeff, class GridFct>
    void operator()(Coeff& coeff, GridFct const& gf) const
    {
      (*this)(coeff, gf, FakeContainer<bool,true>{});
    }
 
    SimpleInterpolator(Basis const& basis, TreePath treePath = {})
      : basis_{basis}
      , treePath_{treePath}
    {}
 
    Basis const& basis_;
    TreePath treePath_ = {};
  };
 
 
  template <class Tag>
  struct InterpolatorFactory;
 
  template <>
  struct InterpolatorFactory<tag::assign>
  {
    template <class Basis,
              class TreePath = Dune::TypeTree::TreePath<>>
    static auto create(Basis const& basis, TreePath treePath = {})
    {
      return SimpleInterpolator<Basis,TreePath>{basis, treePath};
    }
  };
 
} // end namespace AMDiS