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>

namespace AMDiS
{
  namespace tag
  {
    struct assign {};
  }

  template <class Basis,
            class TreePath = Dune::TypeTree::HybridTreePath<>>
  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::HybridTreePath<>>
    static auto create(Basis const& basis, TreePath treePath = {})
    {
      return SimpleInterpolator<Basis,TreePath>{basis, treePath};
    }
  };

} // end namespace AMDiS