VariadicGridFunctionOperator.hpp

#pragma once
 
#include <cassert>
#include <tuple>
#include <type_traits>
 
#include <dune/common/hybridutilities.hh>
 
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/Order.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
#include <amdis/GridFunctionOperatorTransposed.hpp>
#include <amdis/GridFunctionOperator.hpp>
#include <amdis/Output.hpp>
#include <amdis/utility/QuadratureFactory.hpp>
 
namespace AMDiS
{
  template <class Imp, class... LF>
  class VariadicGridFunctionLocalOperator;
 
 
 
  template <class Imp, class... GF>
  class VariadicGridFunctionOperator
  {
  public:
    using Implementation = Imp;
 
 
    template <class... GridFct, class Impl>
    VariadicGridFunctionOperator(Impl&& impl, int derivDeg,
                                 int gridFctOrder, GridFct&&... gridFct)
      : gridFct_(std::forward_as_tuple(gridFct...))
      , impl_(FWD(impl))
      , derivDeg_(derivDeg)
      , gridFctOrder_(gridFctOrder)
    {}
 
    template <class GridView>
    void update(GridView const&) { /* do nothing */ }
 
    friend auto localOperator(VariadicGridFunctionOperator const& op)
    {
      return std::apply([&op](auto const& ... gfs){
        return VariadicGridFunctionLocalOperator{op.impl_,
          op.derivDeg_, op.gridFctOrder_, localFunction(gfs)...};
      }, op.gridFct_);
    }
 
  private:
    std::tuple<GF...> gridFct_;
 
    Implementation impl_;
 
    int derivDeg_;
 
    int gridFctOrder_;
  };
 
  template <class... GridFct, class Impl>
  VariadicGridFunctionOperator(Impl const& impl, int, int, GridFct const&... gridFct)
    -> VariadicGridFunctionOperator< Impl, GridFct... >;
 
 
 
  template <class Imp, class... LF>
  class VariadicGridFunctionLocalOperator
  {
  private:
 
    using Implementation = Imp;
 
  public:
 
    template <class... LocalFct, class Impl>
    VariadicGridFunctionLocalOperator(Impl&& impl, int derivDeg,
                       int localFctOrder, LocalFct&&... localFct)
      : localFct_(std::forward_as_tuple(localFct...))
      , impl_(FWD(impl))
      , derivDeg_(derivDeg)
      , localFctOrder_(localFctOrder)
    {}
 
 
    template <class Element>
    void bind(Element const& element)
    {
      std::apply([&element](auto&... lf){
        (lf.bind(element),...);
      },
      localFct_);
    }
 
    void unbind()
    {
      std::apply([](auto&... lf){
        (lf.unbind(),...);
      },
      localFct_);    }
 
 
    template <class CG, class RN, class CN, class Mat>
    void assemble(CG const& contextGeo, RN const& rowNode, CN const& colNode,
                  Mat& elementMatrix) const
    {
      auto const& quad = getQuadratureRule(contextGeo.geometry(),
        derivDeg_, localFctOrder(), rowNode, colNode);
 
        std::apply([&](auto const& ... lf){
          impl().assemble(contextGeo, rowNode, colNode, quad, lf..., elementMatrix);
        }
        ,localFct_);
    }
 
 
    template <class CG, class Node, class Vec>
    void assemble(CG const& contextGeo, Node const& node,
                  Vec& elementVector) const
    {
      auto const& quad = getQuadratureRule(contextGeo.geometry(),
        derivDeg_, localFctOrder(), node);
 
      std::apply([&](auto const& ... lf){
        impl().assemble(contextGeo, node, quad, lf..., elementVector);
      }, localFct_);
    }
 
    Implementation      & impl()       { return impl_; }
    Implementation const& impl() const { return impl_; }
 
  protected:
    // Return the order of the coefficients function, either localFctOrder_ or
    // the order of the local-function
    int localFctOrder() const
    {
      if (localFctOrder_ >= 0)
        return localFctOrder_;
      else {
        int coeffOrder = -1;
        // Try to get polynomial order of local-functions
        std::apply([&coeffOrder, this](auto const& ... lf){
          coeffOrder = std::max({coeffOrder, getOrder(lf)...});
        }, localFct_);
        test_exit(coeffOrder >= 0,
          "Polynomial degree of coefficients cannot be determined. "
          "Please provide a quadrature order manually.");
 
        return coeffOrder;
      }
    }
 
    // get the polynomial order of a single local function or -1 if it cannot be determined.
    template<class LFF>
    constexpr int getOrder(LFF const& lf) const{
      if constexpr (Concepts::Polynomial<LFF>)
        return order(lf);
      else
        return -1;
    }
 
  private:
    std::tuple<LF...> localFct_;
 
    Implementation impl_;
 
    int derivDeg_;
 
    int localFctOrder_;
  };
 
  // deduction guide
  template <class Impl, class... LocalFct>
  VariadicGridFunctionLocalOperator(Impl const& impl, int, int, LocalFct const&... localFct)
    -> VariadicGridFunctionLocalOperator<Impl, LocalFct...>;
 
 
  template <class Tag, class LocalContext>
  struct VariadicGridFunctionOperatorRegistry {};
 
  template <class Tag, class... Expr>
  struct VariadicOperatorTerm
  {
    Tag tag;
    std::tuple<Expr...> expr;
    int gridFctDeg;
 
    VariadicOperatorTerm(Tag tag,int gridFctDeg, Expr const&... expr)
      : tag(tag)
      , expr(std::make_tuple(expr...))
      , gridFctDeg(gridFctDeg)
    {}
  };
 
 
  template <class Tag, class... Expr>
  auto makeOperator(Tag const& tag, int gridFctDeg, Expr const&... expr)
        ->std::enable_if_t<sizeof...(Expr) != 1, VariadicOperatorTerm<Tag, Expr...>>
  {
    return VariadicOperatorTerm{tag, gridFctDeg, expr...};
  }
 
 
#ifndef DOXYGEN
 
 
  template <class Context, class Tag, class... Expr, class GridView>
  auto makeOperator(VariadicOperatorTerm<Tag,Expr...> const& op, GridView const& gridView)
  {
    using Registry = VariadicGridFunctionOperatorRegistry<Tag,Context>;
    if constexpr (Dune::Std::is_detected_v<IsTransposed,Registry,Tag>) {
      auto impl = typename Registry::type{Registry::transposedTag(op.tag)};
      return GridFunctionOperatorTransposed{
        std::apply([&](auto const&... expr){
          return VariadicGridFunctionOperator{
                  std::move(impl),
                  Registry::degree,
                  op.gridFctDeg,
                  makeGridFunction(expr,gridView)...} ;
        }, op.expr)
      };
    } else {
      auto impl = typename Registry::type{op.tag};
      return VariadicGridFunctionOperator{
        std::apply([&](auto const&... expr){
          return VariadicGridFunctionOperator{
                  std::move(impl),
                  Registry::degree,
                  op.gridFctDeg,
                  makeGridFunction(expr,gridView)...} ;
        }, op.expr)
      };
    }
  }
 
#endif // DOXYGEN
 
} // end namespace AMDiS