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