AnalyticGridFunction.hpp

#pragma once

#include <optional>
#include <type_traits>

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

#include <amdis/Operations.hpp>
#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/Order.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>

namespace AMDiS
{
#ifndef DOXYGEN
  template <class Signature, class LocalContext, class Function>
  class AnalyticLocalFunction;
#endif

  template <class R, class D, class LC, class Function>
  class AnalyticLocalFunction<R(D), LC, Function>
  {
  public:
    using Domain = D;

    using Range = R;

    using LocalContext = LC;

    enum { hasDerivative = true };

  private:
    using Geometry = typename LC::Geometry;

  public:
    AnalyticLocalFunction(Function const& fct)
      : fct_{fct}
    {}

    void bind(LocalContext const& element)
    {
      localContext_.emplace(element);
      geometry_.emplace(element.geometry());
    }

    void unbind()
    {
      geometry_.reset();
      localContext_.reset();
    }
    bool bound() const 
    {
      return !!localContext_;
    }

    LocalContext const& localContext() const
    {
      assert( !!localContext_ );
      return *localContext_;
    }

    Range operator()(Domain const& local) const
    {
      assert( !!geometry_ );
      return fct_(geometry_.value().global(local));
    }

    Function const& fct() const
    {
      return fct_;
    }

  private:
    Function fct_;
    std::optional<LocalContext> localContext_;
    std::optional<Geometry> geometry_;
  };



  template <class Sig, class LC, class F,
    REQUIRES(Concepts::HasFunctorOrder<F,1>)>
  int order(AnalyticLocalFunction<Sig,LC,F> const& lf)
  {
    return order(lf.fct(),1);
  }



  template <class R, class D, class LC, class F, class Type>
  auto derivativeOf(AnalyticLocalFunction<R(D),LC,F> const& lf, Type const& type)
  {
    static_assert(Concepts::HasDerivative<F,Type>,
      "No derivative(F,DerivativeType) defined for Functor F of AnalyticLocalFunction.");

    auto df = derivativeOf(lf.fct(), type);

    using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
    using DerivativeSignature = typename DerivativeTraits<RawSignature,Type>::Range(D);
    auto ldf = AnalyticLocalFunction<DerivativeSignature, LC, decltype(df)>{df};
    // bind if lf is bound
    if (lf.bound())
      ldf.bind(lf.localContext());
    return ldf;
  }


  template <class Function, class GridView>
  class AnalyticGridFunction
  {
  public:
    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
    using Domain = typename EntitySet::GlobalCoordinate;
    using Range = remove_cvref_t<std::result_of_t<Function(Domain)>>;

    enum { hasDerivative = true };

  private:
    using Element = typename EntitySet::Element;
    using LocalDomain = typename EntitySet::LocalCoordinate;
    using LocalFunction = AnalyticLocalFunction<Range(LocalDomain), Element, Function>;

  public:
    AnalyticGridFunction(Function const& fct, GridView const& gridView)
      : fct_{fct}
      , entitySet_{gridView}
    {}

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

    friend LocalFunction localFunction(AnalyticGridFunction const& gf)
    {
      return {gf.fct_};
    }

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

    Function const& fct() const { return fct_; }

  private:
    Function fct_;
    EntitySet entitySet_;
  };



  template <class F, class GV, class Type>
  auto derivativeOf(AnalyticGridFunction<F,GV> const& gf, Type const& type)
  {
    static_assert(Concepts::HasDerivative<F,Type>,
      "No derivative(F,DerivativeType) defined for Functor of AnalyticLocalFunction.");

    auto df = derivativeOf(gf.fct(), type);
    return AnalyticGridFunction<decltype(df), GV>{df, gf.entitySet().gridView()};
  }


#ifndef DOXYGEN
  // A pre-GridFunction that just stores the function
  template <class Function>
  struct AnalyticPreGridFunction
  {
    using Self = AnalyticPreGridFunction;

    struct Creator
    {
      template <class GridView>
      static auto create(Self const& self, GridView const& gridView)
      {
        using Coordinate = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
        static_assert(std::is_invocable_v<Function, Coordinate>);

        return AnalyticGridFunction<Function, GridView>{self.fct_, gridView};
      }
    };

    Function fct_;
  };

  namespace Traits
  {
    template <class Functor>
    struct IsPreGridFunction<AnalyticPreGridFunction<Functor>>
      : std::true_type {};
  }
#endif



  template <class Function>
  auto evalAtQP(Function const& f)
  {
    return AnalyticPreGridFunction<Function>{f};
  }

} // end namespace AMDiS