GridFunctionOperator.hpp

#pragma once

#include <cassert>
#include <type_traits>

#include <amdis/GridFunctions.hpp>
#include <amdis/LocalOperator.hpp>
#include <amdis/common/Transposed.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/typetree/FiniteElementType.hpp>
#include <amdis/utility/QuadratureFactory.hpp>

namespace AMDiS
{

  template <class Derived, class LC, class GF>
  class GridFunctionOperatorBase
      : public LocalOperator<Derived, LC>
  {
    template <class, class>
    friend class LocalOperator;

    using ContextType = Impl::ContextTypes<LC>;
    using Super = LocalOperator<Derived, LC>;

  private:
    using GridFunction = GF;

    using LocalFunction = decltype(localFunction(std::declval<GF>()));

    using Element = typename ContextType::Entity;

    using Geometry = typename Element::Geometry;

    using LocalCoordinate = typename GF::EntitySet::LocalCoordinate;

    using QuadFactory = QuadratureFactory<typename Geometry::ctype, LC::mydimension, LocalFunction>;

  public:

    template <class GridFct>
    GridFunctionOperatorBase(GridFct&& gridFct, int termOrder)
      : gridFct_(FWD(gridFct))
      , termOrder_(termOrder)
    {}

    template <class PQF>
    void setQuadFactory(PQF&& pre)
    {
      using ctype = typename Geometry::ctype;
      quadFactory_ = makeUniquePtr(
        makeQuadratureFactory<ctype, LC::mydimension, LocalFunction>(FWD(pre)));
    }

  protected:
    auto coefficient(LocalCoordinate const& local) const
    {
      assert( this->bound_ );
      return (*localFct_)(local);
    }

    template <class... Nodes>
    auto const& getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
    {
      assert( bool(quadFactory_) );
      int quadDegree = this->getDegree(termOrder_, quadFactory_->order(), nodes...);
      return quadFactory_->rule(type, quadDegree);
    }

  private:

    void bind_impl(Element const& element, Geometry const& geometry)
    {
      assert( bool(quadFactory_) );
      localFct_.emplace(localFunction(gridFct_));
      localFct_->bind(element);
      quadFactory_->bind(localFct_.value());
    }

    void unbind_impl()
    {
      localFct_->unbind();
      localFct_.reset();
    }

  private:
    GridFunction gridFct_;

    std::optional<LocalFunction> localFct_;

    std::shared_ptr<QuadFactory> quadFactory_;

    int termOrder_ = 0;
  };


  template <class Derived, class Transposed>
  class GridFunctionOperatorTransposed
      : public LocalOperator<Derived, typename Transposed::LocalContext>
  {
    template <class, class>
    friend class LocalOperator;

    template <class T, class... Args>
    using Constructable = decltype( new T(std::declval<Args>()...) );

  public:
    template <class... Args,
      std::enable_if_t<Dune::Std::is_detected<Constructable, Transposed, Args...>::value, int> = 0>
    GridFunctionOperatorTransposed(Args&&... args)
      : transposedOp_(FWD(args)...)
    {}

    template <class PQF>
    void setQuadFactory(PQF&& pre)
    {
      transposedOp_.setQuadFactory(FWD(pre));
    }

  private:
    template <class Element, class Geometry>
    void bind_impl(Element const& element, Geometry const& geometry)
    {
      transposedOp_.bind(element, geometry);
    }

    void unbind_impl()
    {
      transposedOp_.unbind();
    }


    template <class CG, class RN, class CN, class Mat>
    void getElementMatrix(CG const& contextGeometry, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
    {
      auto elementMatrixTransposed = transposed(elementMatrix);
      transposedOp_.getElementMatrix(
        contextGeometry, colNode, rowNode, elementMatrixTransposed);
    }

  private:
    Transposed transposedOp_;
  };


#ifndef DOXYGEN
  template <class Tag, class PreGridFct, class PQF>
  struct PreGridFunctionOperator
  {
    Tag tag;
    PreGridFct expr;
    PQF quadFactory;
  };
#endif

  template <class Tag, class Expr, class... QuadratureArgs>
  auto makeOperator(Tag tag, Expr&& expr, QuadratureArgs&&... args)
  {
    auto pqf = makePreQuadratureFactory(FWD(args)...);
    return PreGridFunctionOperator<Tag, TYPEOF(expr), TYPEOF(pqf)>{tag, FWD(expr), std::move(pqf)};
  }

#ifndef DOXYGEN


  template <class Tag, class LC, class GridFct>
  class GridFunctionOperator
      : public GridFunctionOperatorBase<GridFunctionOperator<Tag, LC, GridFct>, LC, GridFct>
  {};

  template <class Context, class Tag, class GF, class QF>
  auto makeGridFunctionOperator(Tag tag, GF&& gf, QF&& qf)
  {
    GridFunctionOperator<Tag, Context, TYPEOF(gf)> gfo{tag, FWD(gf)};
    gfo.setQuadFactory(FWD(qf));
    return gfo;
  }


  template <class Context, class... Args, class GridView>
  auto makeLocalOperator(PreGridFunctionOperator<Args...> op, GridView const& gridView)
  {
    auto gf = makeGridFunction(std::move(op.expr), gridView);
    return makeGridFunctionOperator<Context>(op.tag, std::move(gf), std::move(op.quadFactory));
  }

  template <class Context, class... Args, class GridView>
  auto makeLocalOperator(std::reference_wrapper<PreGridFunctionOperator<Args...>> op_ref, GridView const& gridView)
  {
    PreGridFunctionOperator<Args...> const& op = op_ref;
    auto gf = makeGridFunction(std::ref(op.expr), gridView);
    return makeGridFunctionOperator<Context>(op.tag, std::move(gf), op.quadFactory);
  }

#endif // DOXYGEN

} // end namespace AMDiS