QuadratureFactory.hpp

#pragma once

#include <dune/common/exceptions.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>

#include <amdis/Output.hpp>
#include <amdis/common/Order.hpp>
#include <amdis/functions/Order.hpp>

namespace AMDiS
{
  template <class Geometry, class RN, class CN>
  int getQuadratureDegree(Geometry const& geometry, int derivOrder, int coeffDegree,
                          RN const& rowNode, CN const& colNode)
  {
    test_exit(coeffDegree >= 0,
      "Polynomial degree of coefficients cannot be determined. "
      "Please provide a quadrature order manually.");

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

    int degree = psiDegree + phiDegree + coeffDegree;
    if (geometry.type().isSimplex())
      degree -= derivOrder;
    if (geometry.affine())
      degree += 1;

    return degree;
  }

  template <class Geometry, class Node>
  int getQuadratureDegree(Geometry const& geometry, int derivOrder, int coeffDegree,
                          Node const& node)
  {
    test_exit(coeffDegree >= 0,
      "Polynomial degree of coefficients cannot be determined. "
      "Please provide a quadrature order manually.");

    int psiDegree = order(node);
    int degree = psiDegree + coeffDegree;

    if (geometry.type().isSimplex())
      degree -= derivOrder;
    if (geometry.affine())
      degree += 1;

    return degree;
  }

  template <class Geometry, class... Nodes>
  auto const& getQuadratureRule(Geometry const& geometry, int derivOrder,
                                int coeffDegree, Nodes const&... nodes)
  {
    int degree = getQuadratureDegree(geometry, derivOrder, coeffDegree, nodes...);
    using QuadRules
      = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>;
    return QuadRules::rule(geometry.type(), degree);
  }

} // end namespace AMDiS