[!127] Added a new quadraturetype.

Merge branch 'master' into 'master'

ref:core/dune-geometry GaussJacobiArbitraryOrder is a quadrature type that
respects the weight function resulting from conical product.

The rules are computed on demand and currently have a "maximal" order of 127.
This upper limit is arbitrary and can be changed.

See merge request [core/dune-geometry!127]

  [core/dune-geometry!127]: gitlab.dune-project.org/core/dune-geometry/merge_requests/127

CHANGELOG.md

@@ -10,6 +10,15 @@
   `GeometryTypes::cube(dim)` to construct one- or two-dimensional `GeometryType` objects.
 - Geometry implementations now export a type `Volume` that is used for the return
   value of the `volume` methods.  So does the generic `ReferenceElement` implementation.
+-   More efficient quadrature rules for simplices are available that
+    need less quadrature points to achive the same order.  For now these
+    have to be explicitly requested:
+    ```c++
+    auto&& rule = Dune::QuadratureRules<...>::rule(..., Dune::QuadratureType::GaussJacobi_n_0);
+    ```
+    See [!127].
+
+    [!127]: https://gitlab.dune-project.org/core/dune-geometry/merge_requests/127
 
 # Release 2.6
 

dune/geometry/quadraturerules.hh

@@ -78,11 +78,41 @@ namespace Dune {
    */
   namespace QuadratureType {
     enum Enum {
+      /** \brief Gauss-Legendre rules (default)
+      *
+      *  -1D: Gauss-Jacobi rule with parameters \f$\alpha = \beta =0 \f$
+      *  -higher dimension: For the 2D/3D case efficient rules for certain geometries may be used if available.
+      *                     Higher dimensional quadrature rules are constructed via \p TensorProductQuadratureRule.
+      *                     In this case the 1D rules eventually need higher order to compensate occuring weight functions(i.e. simplices).
+      */
       GaussLegendre = 0,
 
+      /** \brief Gauss-Jacobi rules with \f$\alpha =1\f$
+      *
+      *  -1D Gauss-Jacobi rule with parameters \f$\alpha =1,\ \beta =0 \f$
+      *  -Is used to construct efficient simplex quadrature rules of higher order
+      */
       GaussJacobi_1_0 = 1,
+
+      /** \brief Gauss-Legendre rules with \f$\alpha =2\f$
+      *
+      *  -1D Gauss-Jacobi rule with parameters \f$\alpha =2,\ \beta =0 \f$
+      *  -Is used to construct efficient simplex quadrature rules of higher order
+      */
       GaussJacobi_2_0 = 2,
 
+      /** \brief Gauss-Legendre rules with \f$\alpha =n\f$.
+      *
+      *  -1D: Gauss-Jacobi rule with parameters \f$\alpha = n,\ \beta =0 \f$
+      *  -higher dimension: For the 2D/3D case efficient rules for certain geometries may be used if available.
+      *                     Higher dimensional quadrature rules are constructed via \p TensorProductQuadratureRule.
+      *                     In this case the 1D rules respect eventually occuring weight functions(i.e. simplices).
+      *  -The rules for high dimension or order are computed at run time and only floating point number types are supported.(LAPACK is needed for this case)
+      *  -Most efficient quadrature type for simplices.
+      *
+      *   \note For details please use the book "Approximate Calculation of Multiple Integrals" by A.H. Stroud published in 1971.
+      */
+      GaussJacobi_n_0 = 3,
       GaussLobatto = 4,
       size
     };
@@ -741,6 +771,8 @@ namespace Dune {
           return Jacobi2QuadratureRule1D<ctype>::highest_order;
         case QuadratureType::GaussLobatto :
           return GaussLobattoQuadratureRule1D<ctype>::highest_order;
+        case QuadratureType::GaussJacobi_n_0 :
+          return JacobiNQuadratureRule1D<ctype>::maxOrder();
         default :
           DUNE_THROW(Exception, "Unknown QuadratureType");
         }
@@ -760,6 +792,8 @@ namespace Dune {
           return Jacobi2QuadratureRule1D<ctype>(p);
         case QuadratureType::GaussLobatto :
           return GaussLobattoQuadratureRule1D<ctype>(p);
+        case QuadratureType::GaussJacobi_n_0 :
+          return JacobiNQuadratureRule1D<ctype>(p);
         default :
           DUNE_THROW(Exception, "Unknown QuadratureType");
         }
@@ -785,7 +819,7 @@ namespace Dune {
     static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
     {
       if (t.isSimplex()
-        && qt == QuadratureType::GaussLegendre
+        && ( qt == QuadratureType::GaussLegendre || qt == QuadratureType::GaussJacobi_n_0 )
         && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
       {
         return SimplexQuadratureRule<ctype,dim>(p);
@@ -813,8 +847,9 @@ namespace Dune {
     }
     static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
     {
+
       if (t.isSimplex()
-        && qt == QuadratureType::GaussLegendre
+        && ( qt == QuadratureType::GaussLegendre || qt == QuadratureType::GaussJacobi_n_0 )
         && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
       {
         return SimplexQuadratureRule<ctype,dim>(p);

dune/geometry/quadraturerules/CMakeLists.txt

@@ -7,6 +7,7 @@ install(FILES
   gausslobatto_imp.hh
   jacobi_1_0_imp.hh
   jacobi_2_0_imp.hh
+  jacobi_n_0.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/geometry/quadraturerules)
 
 exclude_from_headercheck(

dune/geometry/quadraturerules/jacobi_n_0.hh

+#ifndef DUNE_GEOMETRY_QUADRATURERULES_JACOBI_N_0_H
+#define DUNE_GEOMETRY_QUADRATURERULES_JACOBI_N_0_H
+
+#include <vector>
+#include <type_traits>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/exceptions.hh>
+
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+#include <dune/common/dynmatrixev.hh>
+
+
+namespace Dune {
+
+  /*
+   * This class calculates 1D quadrature rules of dynamic order, that respect the
+   * weightfunctions resulting from the conical product.
+   * For quadrature rules with order <= highest order (61) and dimension <= 3 exact quadrature rules are used.
+   * Else the quadrature rules are computed by Lapack and thus a floating point type is required.
+   * For more information see the comment in `tensorproductquadrature.hh`
+   *
+   */
+
+  template<typename ct>
+  class JacobiNQuadratureRule1D : public QuadratureRule<ct,1>
+  {
+  public:
+    // compile time parameters
+    enum { dim=1 };
+
+  private:
+
+    typedef QuadratureRule<ct, dim> Rule;
+
+    friend class QuadratureRuleFactory<ct,dim>;
+
+    template< class ctype, int dimension>
+    friend class TensorProductQuadratureRule;
+
+    explicit JacobiNQuadratureRule1D (int const order, int const alpha=0)
+      : Rule( GeometryTypes::line )
+    {
+      if (unsigned(order) > maxOrder())
+        DUNE_THROW(QuadratureOrderOutOfRange, "Quadrature rule " << order << " not supported!");
+      auto&& rule = decideRule(order,alpha);
+      for( auto qpoint : rule )
+        this->push_back(qpoint);
+      this->delivered_order = 2*rule.size()-1;
+
+    }
+
+    static unsigned maxOrder()
+    {
+      return 127; // can be changed
+    }
+
+    QuadratureRule<ct,1> decideRule(int const degree, int const alpha)
+    {
+      const auto maxOrder = std::min(   unsigned(GaussQuadratureRule1D<ct>::highest_order),
+                              std::min( unsigned(Jacobi1QuadratureRule1D<ct>::highest_order),
+                                        unsigned(Jacobi2QuadratureRule1D<ct>::highest_order))
+                              );
+      return unsigned(degree) < maxOrder ? decideRuleExponent(degree,alpha) : UseLapackOrError<ct>(degree, alpha);
+    }
+
+    template<typename type, std::enable_if_t<!(std::is_floating_point<type>::value)>* = nullptr>
+    QuadratureRule<ct,1> UseLapackOrError( int const degree, int const alpha)
+    {
+      DUNE_THROW(QuadratureOrderOutOfRange, "Quadrature rule with degree: "<< degree << " and jacobi exponent: "<< alpha<< " is not supported for this type!");
+    }
+
+    template<typename type, std::enable_if_t<std::is_floating_point<type>::value>* = nullptr>
+    QuadratureRule<ct,1> UseLapackOrError( int const degree, int const alpha)
+    {
+      return jacobiNodesWeights<ct>(degree, alpha);
+    }
+
+    QuadratureRule<ct,1> decideRuleExponent(int const degree, int const alpha)
+    {
+      switch(alpha)
+      {
+        case 0 :  return QuadratureRules<ct,1>::rule(GeometryTypes::line, degree, QuadratureType::GaussLegendre);
+        case 1 :  {
+                    // scale already existing weights by 0.5
+                    auto&& rule = QuadratureRules<ct,1>::rule(GeometryTypes::line, degree, QuadratureType::GaussJacobi_1_0);
+                    QuadratureRule<ct,1> quadratureRule;
+                    quadratureRule.reserve(rule.size());
+                    for( auto qpoint : rule )
+                      quadratureRule.push_back(QuadraturePoint<ct,1>(qpoint.position(),ct(0.5)*qpoint.weight()));
+                    return quadratureRule;
+                  }
+        case 2 :  {
+                    // scale already existing weights by 0.25
+                    auto&& rule = QuadratureRules<ct,1>::rule(GeometryTypes::line, degree, QuadratureType::GaussJacobi_2_0);
+                    QuadratureRule<ct,1> quadratureRule;
+                    quadratureRule.reserve(rule.size());
+                    for( auto qpoint : rule )
+                      quadratureRule.push_back(QuadraturePoint<ct,1>(qpoint.position(),ct(0.25)*qpoint.weight()));
+                    return quadratureRule;
+                  }
+        default : return UseLapackOrError<ct>(degree,alpha);
+      }
+    }
+
+    // computes the nodes and weights for the weight function (1-x)^alpha, which are exact for given polynomials with degree "degree"
+    template<typename ctype, std::enable_if_t<std::is_floating_point<ctype>::value>* = nullptr>
+    QuadratureRule<ct,1> jacobiNodesWeights(int const degree, int const alpha)
+    {
+      using std::sqrt;
+
+      typedef DynamicVector<ct> Vector;
+      typedef DynamicMatrix<ct> Matrix;
+
+      // compute the degree of the needed jacobi polymomial
+      const int n = degree/2 +1;
+
+      Matrix J(n,n,0);
+
+      J[0][0] = (ct) -alpha/(2 + alpha);
+      for(int i=1; i<n; ++i)
+      {
+        ct a_i = (ct) (2*i*(i + alpha)) /( (2*i + alpha - 1)*(2*i + alpha) );
+        ct c_i = (ct) (2*i*(i + alpha)) /( (2*i + alpha + 1)*(2*i + alpha) );
+
+        J[i][i] =  (ct) - alpha*alpha /( (2*i + alpha + 2)*(2*i + alpha) );
+        J[i][i-1] = sqrt(a_i*c_i);
+        J[i-1][i] = J[i][i-1];
+      }
+
+      DynamicVector<std::complex<double> > eigenValues(n,0);
+      std::vector<DynamicVector<double> > eigenVectors(n, DynamicVector<double>(n,0));
+
+      DynamicMatrixHelp::eigenValuesNonSym(J, eigenValues, &eigenVectors);
+
+      ct mu = (ct) 1/(alpha + 1);
+
+      QuadratureRule<ct,1> quadratureRule;
+      quadratureRule.reserve(n);
+      for (int i=0; i<n; ++i)
+      {
+        auto&& eV0 = eigenVectors[i][0];
+        ct weight =  mu * eV0*eV0;
+        Vector node(1,0.5*eigenValues[i].real() + 0.5);
+
+        // bundle the nodes and the weights
+        QuadraturePoint<ct,1> temp(node, weight);
+        quadratureRule.push_back(temp);
+      }
+
+      return quadratureRule;
+
+    }
+
+  };
+
+}
+
+#endif

dune/geometry/quadraturerules/tensorproductquadrature.hh

@@ -7,6 +7,7 @@
 #include <bitset>
 
 #include <dune/geometry/type.hh>
+#include <dune/geometry/quadraturerules/jacobi_n_0.hh>
 
 namespace Dune
 {
@@ -108,8 +109,16 @@ namespace Dune
     void conicalProduct(const BaseQuadrature & baseQuad, unsigned int order, QuadratureType::Enum qt)
     {
       typedef QuadratureRule<ctype,1> OneDQuadrature;
-      const OneDQuadrature & onedQuad =
-        QuadratureRules<ctype,1>::rule(GeometryTypes::line, order + dim-1, qt);
+
+      OneDQuadrature onedQuad;
+      bool respectWeightFunction = false;
+      if( qt != QuadratureType::GaussJacobi_n_0)
+        onedQuad = QuadratureRules<ctype,1>::rule(GeometryTypes::line, order + dim-1, qt);
+      else
+      {
+        onedQuad = JacobiNQuadratureRule1D<ctype>(order,dim-1);
+        respectWeightFunction = true;
+      }
 
       const unsigned int baseQuadSize = baseQuad.size();
       for( unsigned int bqi = 0; bqi < baseQuadSize; ++bqi )
@@ -128,8 +137,12 @@ namespace Dune
             point[ i ] = scale * basePoint[ i ];
 
           ctype weight = baseWeight * onedQuad[oqi].weight( );
-          for ( unsigned int p = 0; p<dim-1; ++p)
-            weight *= scale;                    // pow( scale, dim-1 );
+          if (!respectWeightFunction)
+          {
+            for ( unsigned int p = 0; p<dim-1; ++p)
+              weight *= scale;                    // pow( scale, dim-1 );
+          }
+
           this->push_back( QPoint(point, weight) );
         }
       }
@@ -146,9 +159,12 @@ namespace Dune
       if (isPrism)
         order = std::min
           (order, QuadratureRules<ctype,1>::maxOrder(GeometryTypes::line, qt));
-      else
+      else if (qt != QuadratureType::GaussJacobi_n_0)
         order = std::min
           (order, QuadratureRules<ctype,1>::maxOrder(GeometryTypes::line, qt)-(dim-1));
+      else
+        order = std::min
+          (order, QuadratureRules<ctype,1>::maxOrder(GeometryTypes::line, qt));
       return order;
     }
 

dune/geometry/test/test-quadrature.cc

@@ -220,6 +220,7 @@ int main (int argc, char** argv)
     check<double,4>(Dune::GeometryTypes::cube(4), std::min(maxOrder, unsigned(31)),
                     Dune::QuadratureType::GaussLobatto);
     check<double,4>(Dune::GeometryTypes::simplex(4), maxOrder);
+    check<double,4>(Dune::GeometryTypes::simplex(4), maxOrder, Dune::QuadratureType::GaussJacobi_n_0);
     check<double,3>(Dune::GeometryTypes::prism, maxOrder);
     check<double,3>(Dune::GeometryTypes::pyramid, maxOrder);