LocalOperator.hpp

#pragma once

#include <cassert>
#include <type_traits>

// #include <amdis/GridFunctions.hpp>
#include <amdis/ContextGeometry.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/functions/Order.hpp>
#include <amdis/typetree/FiniteElementType.hpp>

namespace AMDiS
{

  template <class Derived, class LC>
  class LocalOperator
  {
    using ContextType = Impl::ContextTypes<LC>;

  public:
    using LocalContext = LC;

    using Element = typename ContextType::Entity;

    using Geometry = typename Element::Geometry;

    template <class GridView>
    void init(GridView const& gridView)
    {
      derived().init_impl(gridView);
    }


    void bind(Element const& element, Geometry const& geometry)
    {
      if (!bound_) {
        isSimplex_ = geometry.type().isSimplex();
        isAffine_ = geometry.affine();
        bound_ = true;
      }

      derived().bind_impl(element, geometry);
    }

    void unbind()
    {
      derived().unbind_impl();
      bound_ = false;
    }


    template <class CG, class RN, class CN, class Mat>
    void calculateElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
    {
      assert( bound_ );
      derived().getElementMatrix(contextGeo, rowNode, colNode, elementMatrix);
    }


    template <class CG, class Node, class Vec>
    void calculateElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
    {
      assert( bound_ );
      derived().getElementVector(contextGeo, node, elementVector);
    }


    Derived      & derived()       { return static_cast<Derived&>(*this); }
    Derived const& derived() const { return static_cast<Derived const&>(*this); }


  protected:

    // default implementation. Can be overridden in the derived classes
    template <class GridView>
    void init_impl(GridView const& /*gridView*/) {}

    // default implementation. Can be overridden in the derived classes
    template <class Element, class Geometry>
    void bind_impl(Element const& /*element*/, Geometry const& /*geometry*/) {}

    // default implementation. Can be overridden in the derived classes
    void unbind_impl() {}


    template <class RN, class CN>
    int getDegree(int derivOrder, int coeffDegree, RN const& rowNode, CN const& colNode) const
    {
      assert( bound_ );
      test_warning(coeffDegree >= 0,
        "polynomial order of coefficient function not determined. Use order 4 by default.");

      int psiDegree = order(rowNode);
      int phiDegree = order(colNode);

      int degree = psiDegree + phiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
      if (isSimplex_)
        degree -= derivOrder;
      if (isAffine_)
        degree += 1; // oder += (derivOrder+1)

      return degree;
    }


    template <class Node>
    int getDegree(int derivOrder, int coeffDegree, Node const& node) const
    {
      assert( bound_ );
      test_warning(coeffDegree >= 0,
        "polynomial order of coefficient function not determined. Use order 4 by default.");

      int psiDegree = order(node);

      int degree = psiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
      if (isSimplex_)
        degree -= derivOrder;
      if (isAffine_)
        degree += 1; // oder += (derivOrder+1)

      return degree;
    }

  protected:

    bool isSimplex_ = false;    //< the bound element is a simplex
    bool isAffine_ = false;     //< the bound geometry is affine
    bool bound_ = false;        //< an element is bound to the operator
  };



  template <class Derived, class LC, class GV>
  auto makeLocalOperator(LocalOperator<Derived, LC> const& localOp, GV const& /*gridView*/)
  {
    return localOp.derived();
  }

} // end namespace AMDiS