GridFunctionOperator.hpp

#pragma once
 
#include <cassert>
#include <type_traits>
 
#include <amdis/GridFunctionOperatorTransposed.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/Order.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/utility/QuadratureFactory.hpp>
 
namespace AMDiS
{
  template <class LF, class Imp>
  class GridFunctionLocalOperator;
 
 
 
  template <class GF, class Imp>
  class GridFunctionOperator
  {
  public:
    using GridFunction = GF;
    using Implementation = Imp;
 
 
    template <class GridFct, class Impl>
    GridFunctionOperator(GridFct&& gridFct, Impl&& impl,
                         int derivDeg, int gridFctOrder)
      : gridFct_(FWD(gridFct))
      , impl_(FWD(impl))
      , derivDeg_(derivDeg)
      , gridFctOrder_(gridFctOrder)
    {}
 
    template <class GridView>
    void update(GridView const&) { /* do nothing */ }
 
    friend auto localOperator(GridFunctionOperator const& op)
    {
      return GridFunctionLocalOperator{localFunction(op.gridFct_), op.impl_,
        op.derivDeg_, op.gridFctOrder_};
    }
 
  private:
    GridFunction gridFct_;
 
    Implementation impl_;
 
    int derivDeg_;
 
    int gridFctOrder_;
  };
 
  template <class GridFct, class Impl>
  GridFunctionOperator(GridFct const& gridFct, Impl const& impl, int, int)
    -> GridFunctionOperator<GridFct, Impl>;
 
 
 
  template <class LF, class Imp>
  class GridFunctionLocalOperator
  {
  private:
    using LocalFunction = LF;
 
    using Implementation = Imp;
 
  public:
 
    template <class LocalFct, class Impl>
    GridFunctionLocalOperator(LocalFct&& localFct, Impl&& impl,
                              int derivDeg, int localFctOrder)
      : localFct_(FWD(localFct))
      , impl_(FWD(impl))
      , derivDeg_(derivDeg)
      , localFctOrder_(localFctOrder)
    {}
 
 
    template <class Element>
    void bind(Element const& element)
    {
      localFct_.bind(element);
    }
 
    void unbind()
    {
      localFct_.unbind();
    }
 
 
    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);
      impl().assemble(contextGeo, rowNode, colNode, quad, localFct_, elementMatrix);
    }
 
 
    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);
      impl().assemble(contextGeo, node, quad, localFct_, elementVector);
    }
 
    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
        if constexpr (Concepts::Polynomial<LF>)
          coeffOrder = order(localFct_);
 
        test_exit(coeffOrder >= 0,
          "Polynomial degree of coefficients cannot be determined. "
          "Please provide a quadrature order manually.");
 
        return coeffOrder;
      }
    }
 
  private:
    LocalFunction localFct_;
 
    Implementation impl_;
 
    int derivDeg_;
 
    int localFctOrder_;
  };
 
  // deduction guide
  template <class LocalFct, class Impl>
  GridFunctionLocalOperator(LocalFct const& localFct, Impl const& impl, int, int)
    -> GridFunctionLocalOperator<LocalFct, Impl>;
 
 
  template <class Tag, class LocalContext>
  struct GridFunctionOperatorRegistry {};
 
  template <class Tag, class Expr>
  struct OperatorTerm
  {
    Tag tag;
    Expr expr;
    int gridFctDeg;
 
    OperatorTerm(Tag tag, Expr const& expr, int gridFctDeg = -1)
      : tag(tag)
      , expr(expr)
      , gridFctDeg(gridFctDeg)
    {}
  };
 
  template <class Tag, class Expr>
  auto makeOperator(Tag const& tag, Expr&& expr, int gridFctDeg = -1)
  {
    return OperatorTerm{tag, FWD(expr), gridFctDeg};
  }
 
  template <class Tag, class Expr>
  auto operatorTerm(Tag const& tag, Expr&& expr, int gridFctDeg = -1)
  {
    return OperatorTerm{tag, FWD(expr), gridFctDeg};
  }
 
#ifndef DOXYGEN
 
  template <class R, class Tag>
  using IsTransposed = decltype(R::transposedTag(std::declval<Tag>()));
 
  template <class Context, class Tag, class Expr, class GridView>
  auto makeOperator(OperatorTerm<Tag,Expr> const& op, GridView const& gridView)
  {
    auto gf = makeGridFunction(std::move(op.expr), gridView);
 
    using Registry = GridFunctionOperatorRegistry<Tag,Context>;
    if constexpr (Dune::Std::is_detected_v<IsTransposed,Registry,Tag>) {
      auto impl = typename Registry::type{Registry::transposedTag(op.tag)};
      return GridFunctionOperatorTransposed{
        GridFunctionOperator{std::move(gf), std::move(impl), Registry::degree, op.gridFctDeg}
      };
    } else {
      auto impl = typename Registry::type{op.tag};
      return GridFunctionOperator{std::move(gf), std::move(impl), Registry::degree, op.gridFctDeg};
    }
  }
 
#endif // DOXYGEN
 
} // end namespace AMDiS