CoarsenedGridFunction.hpp

#pragma once

#include <functional>
#include <type_traits>

#include <dune/common/referencehelper.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>

#include <amdis/common/TypeTraits.hpp>
#include <amdis/Output.hpp>

// backport of dune/functions/gridfunctions/coarsenedgridfunction.hh

namespace AMDiS {


template <class GV, class GF>
class CoarsenedGridFunction
{
public:
  using EntitySet = Dune::Functions::GridViewEntitySet<GV,0>;

  using Domain = typename EntitySet::GlobalCoordinate;

  using Range = std::invoke_result_t<GF, Domain>;

  using LocalDomain = typename EntitySet::LocalCoordinate;

  using Element = typename EntitySet::Element;

  using Grid = typename GV::Grid;
  using Geometry = typename Element::Geometry;

  template <class ES, class LF>
  class CoarsenedLocalFunction
  {
  public:
    using Domain = LocalDomain;
    using Range = CoarsenedGridFunction::Range;
    using Element = CoarsenedGridFunction::Element;

    using Mapping = std::function<Domain(Domain const&)>;
    struct ChildElement {
      typename Element::HierarchicIterator it;
      Mapping local;
    };

  public:
    CoarsenedLocalFunction(ES const& entitySet, LF const& localFct, int maxLevel)
      : entitySet_(entitySet)
      , localFct_(localFct)
      , maxLevel_(maxLevel)
    {}

    CoarsenedLocalFunction(ES const& entitySet, LF&& localFct, int maxLevel)
      : entitySet_(entitySet)
      , localFct_(std::move(localFct))
      , maxLevel_(maxLevel)
    {}


    void bind(Element const& element)
    {
      element_.emplace(element);

      childs_.clear();
      if (entitySet_.contains(*element_)) {
        localFct_.bind(*element_);
        return;
      }

      test_exit_dbg(!element_->isLeaf(), "Element is leaf. Cannot traverse its children.");
      for (auto it = element_->hbegin(maxLevel_); it != element_->hend(maxLevel_); ++it) {
        if (entitySet_.contains(*it))
          childs_.emplace_back(ChildElement{.it=it, .local=makeMapping(*element_,*it)});
      }

      test_exit_dbg(!childs_.empty(), "No child element in entitySet found!");
    }

    void unbind()
    {
      element_.reset();
    }

    bool bound() const
    {
      return !!element_;
    }

    Range operator()(Domain const& x) const
    {
      if (childs_.empty())
        return localFct_(x);

      // calculate checkInsideTolerance
      typename Grid::ctype const checkInsideTolerance = std::sqrt(std::numeric_limits<typename Grid::ctype>::epsilon());
      for (auto const& child : childs_) {
        auto refElem = referenceElement(*child.it);
        auto local = child.local(x);
        auto refTypeId = refElem.type().id();
        bool isInside = Dune::Geo::Impl::checkInside(refTypeId, Geometry::mydimension, local, checkInsideTolerance);
        if (isInside) {
          localFct_.bind(*child.it);
          return localFct_(local);
        }
      }

      error_exit("No child element with x in child found!");
      return Range{};
    }

    Element const& localContext() const
    {
      assert(!!element_);
      return *element_;
    }

    template <class LF_ = LF>
    auto makeDerivative() const
      -> CoarsenedLocalFunction<ES,TYPEOF(derivative(std::declval<LF_ const&>()))>
    {
      CoarsenedLocalFunction<ES, TYPEOF(derivative(std::declval<LF_ const &>()))> df{
          entitySet_, derivative(localFct_)};

      // bind derivative if this is bound
      if (bound())
        df.bind(*element_);
      return df;
    }

  protected:
    // Construct the coordinate mapping
    Mapping makeMapping(Element const& coarse, Element fine) const
    {
      Mapping map = [](Domain const& x) { return x; };
      while (coarse != fine && fine.hasFather()) {
        map = [map, geo=fine.geometryInFather()](Domain const& x) {
          return map(geo.local(x));
        };
        fine = fine.father();
      }
      return map;
    }

  private:
    ES const& entitySet_;
    mutable LF localFct_;
    int maxLevel_;
    std::optional<Element> element_;

    std::vector<ChildElement> childs_;
  };



  CoarsenedGridFunction(GV const& gridView, GF const& gridFct)
    : gridView_{gridView}
    , entitySet_{gridView}
    , gridFct_{gridFct}
  {}

  CoarsenedGridFunction(GV const& gridView, GF&& gridFct)
    : gridView_{gridView}
    , entitySet_{gridView}
    , gridFct_{std::move(gridFct)}
  {}

  Range operator()(Domain const& x) const
  {
    return gridFct_(x);
  }

  template <class GF_ = GF>
  auto makeDerivative() const
    -> CoarsenedGridFunction<GV,
        TYPEOF(derivative(Dune::resolveRef(std::declval<GF_ const&>())))>
  {
    return {gridView_, derivative(Dune::resolveRef(gridFct_))};
  }

  auto makeLocalFunction() const
  {
    using LF = TYPEOF(localFunction(Dune::resolveRef(gridFct_)));
    using LocalFunction = CoarsenedLocalFunction<EntitySet,LF>;
    return LocalFunction{gridFct_.entitySet(),
                         localFunction(Dune::resolveRef(gridFct_)),
                         gridView_.grid().maxLevel()};
  }

  EntitySet const& entitySet () const
  {
    return entitySet_;
  }

private:
  GV gridView_;
  EntitySet entitySet_;
  GF gridFct_;
};

} // end namespace AMDiS