DerivativeGridFunction.hpp

#pragma once
 
#include <type_traits>
 
#include <dune/grid/utility/hierarchicsearch.hh>
 
#include <amdis/common/Concepts.hpp>
#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/gridfunctions/Derivative.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
 
namespace AMDiS
{
  namespace Impl
  {
    template <class LF, class Sig>
    struct CheckFunctorConcept
    {
      static_assert(Concepts::Functor<LF, Sig>, "Derivative of LocalFunction can not be called as a functor.");
    };
 
    template <class Traits>
    struct CheckValidRange
    {
      static_assert(!std::is_same_v<typename Traits::Range, Dune::Functions::InvalidRange>, "Invalid Range.");
    };
  }
 
 
  template <class GridFunction, class Type>
  class DerivativeGridFunction
  {
    using GridFctRange = typename GridFunction::Range;
    using GridFctDomain = typename GridFunction::Domain;
    using RawSignature = typename Dune::Functions::SignatureTraits<GridFctRange(GridFctDomain)>::RawSignature;
 
    using Traits = DerivativeTraits<RawSignature, Type>;
    using LocalFunction = TYPEOF( derivativeOf(localFunction(std::declval<GridFunction>()), std::declval<Type>()) ) ;
 
    using LocalFctRange = typename Traits::Range;
    using LocalFctDomain = typename GridFunction::EntitySet::LocalCoordinate;
 
    CHECK_CONCEPT(Impl::CheckValidRange<Traits>);
    CHECK_CONCEPT(Impl::CheckFunctorConcept<LocalFunction, LocalFctRange(LocalFctDomain)>);
 
    enum { hasDerivative = false };
 
  public:
    using Range = typename Traits::Range;
 
    using Domain = GridFctDomain;
 
    using EntitySet = typename GridFunction::EntitySet;
 
  public:
    DerivativeGridFunction(GridFunction const& gridFct, Type const& type)
      : gridFct_{gridFct}
      , type_{type}
    {}
 
    DerivativeGridFunction(GridFunction&& gridFct, Type const& type)
      : gridFct_{std::move(gridFct)}
      , type_{type}
    {}
 
    Range operator()(Domain const& x) const
    {
      auto gv = entitySet().gridView();
 
      using GridView = decltype(gv);
      using Grid = typename GridView::Grid;
      using IS = typename GridView::IndexSet;
 
      Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
 
      auto element = hsearch.findEntity(x);
      auto geometry = element.geometry();
      auto localFct = makeLocalFunction();
      localFct.bind(element);
      return localFct(geometry.local(x));
    }
 
    LocalFunction makeLocalFunction() const
    {
      return derivativeOf(localFunction(gridFct_), type_);
    }
 
    EntitySet const& entitySet() const
    {
      return gridFct_.entitySet();
    }
 
  private:
    GridFunction gridFct_;
    Type type_;
  };
 
 
 
  template <class GridFct, class Type,
            class LocalFct = decltype( localFunction(std::declval<GridFct>()) ),
    REQUIRES(not GridFct::hasDerivative)>
  auto derivativeOf(GridFct const& gridFct, Type const& type)
  {
    static_assert(Concepts::HasDerivative<LocalFct,Type>,
      "derivativeOf(LocalFunction,type) not defined!");
    return DerivativeGridFunction<GridFct,Type>{gridFct, type};
  }
 
 
#ifndef DOXYGEN
  template <class Expr, class Type>
  struct DerivativePreGridFunction
  {
    using Self = DerivativePreGridFunction;
 
    struct Creator
    {
      template <class GridView>
      static auto create(Self const& self, GridView const& gridView)
      {
        return derivativeOf(makeGridFunction(self.expr_, gridView), self.type_);
      }
    };
 
    Expr expr_;
    Type type_ = {};
  };
 
  namespace Traits
  {
    template <class Expr,class Type>
    struct IsPreGridFunction<DerivativePreGridFunction<Expr,Type>>
      : std::true_type {};
  }
#endif
 
 
 
  template <class Expr>
  auto gradientOf(Expr const& expr)
  {
    return DerivativePreGridFunction<Expr, tag::gradient>{expr};
  }
 
  template <class Expr>
  auto divergenceOf(Expr const& expr)
  {
    return DerivativePreGridFunction<Expr, tag::divergence>{expr};
  }
 
  template <class Expr>
  auto partialDerivativeOf(Expr const& expr, std::size_t i)
  {
    return DerivativePreGridFunction<Expr, tag::partial>{expr, tag::partial{i}};
  }
 
} // end namespace AMDiS