LocalAverageInterpolator.hpp

#pragma once

#include <map>
#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 local_average {};
  }

  template <class Basis,
            class TreePath = Dune::TypeTree::HybridTreePath<>>
  struct LocalAverageInterpolator
  {
    template <class Coeff, class GridFct, class BitVector>
    void operator()(Coeff& coeff, GridFct const& gf, BitVector const& bitVec) const
    {
      // storage for values computed during the interpolation
      std::map<typename Basis::MultiIndex, typename Coeff::value_type> values;
      std::map<typename Basis::MultiIndex, std::uint8_t> counter;

      std::vector<typename Coeff::value_type> localCoeff;

      // Obtain a local view of the gridFunction
      auto lf = localFunction(gf);
      auto localView = basis_.localView();
      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)
        {
          // check whether there is a local contribution
          bool any = false;
          for (std::size_t i = 0; i < node.size(); ++i) {
            auto idx = localView.index(node.localIndex(i));
            if(bitVec[idx]) {
              any = true;
              break;
            }
          }

          if (!any)
            return;

          // compute local interpolation coefficients
          node.finiteElement().localInterpolation().interpolate(HierarchicNodeWrapper{tp,lf}, localCoeff);

          // assign relevant local coefficients to temporary storage
          for (std::size_t i = 0; i < node.size(); ++i) {
            auto idx = localView.index(node.localIndex(i));
            if (bitVec[idx]) {
              counter[idx]++;
              values[idx] += localCoeff[i];
            }
          }
        });

        lf.unbind();
      }

      assert(values.size() == counter.size());

      // assign computed values to target vector
      auto value_it = values.begin();
      auto count_it = counter.begin();
      for (; value_it != values.end(); ++value_it, ++count_it) {
        coeff.set(value_it->first, value_it->second/double(count_it->second));
      }
      coeff.finish();
    }

    template <class Coeff, class GridFct>
    void operator()(Coeff& coeff, GridFct const& gf) const
    {
      (*this)(coeff, gf, FakeContainer<bool,true>{});
    }

    LocalAverageInterpolator(Basis const& basis, TreePath treePath = {})
      : basis_{basis}
      , treePath_{treePath}
    {}

    Basis const& basis_;
    TreePath treePath_ = {};
  };


  template <>
  struct InterpolatorFactory<tag::local_average>
  {
    template <class Basis,
              class TreePath = Dune::TypeTree::HybridTreePath<>>
    static auto create(Basis const& basis, TreePath treePath = {})
    {
      return LocalAverageInterpolator<Basis,TreePath>{basis, treePath};
    }
  };

} // end namespace AMDiS