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