]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add QGaussRadau quadrature + tests 15986/head
authorJan Philipp Thiele <thiele@wias-berlin.de>
Fri, 15 Sep 2023 09:30:17 +0000 (11:30 +0200)
committerJan Philipp Thiele <thiele@wias-berlin.de>
Fri, 15 Sep 2023 09:30:17 +0000 (11:30 +0200)
doc/news/changes/minor/20230915Thiele [new file with mode: 0644]
include/deal.II/base/quadrature_lib.h
source/base/quadrature_lib.cc
tests/base/quadrature_test.cc
tests/base/quadrature_test.output

diff --git a/doc/news/changes/minor/20230915Thiele b/doc/news/changes/minor/20230915Thiele
new file mode 100644 (file)
index 0000000..6160f9b
--- /dev/null
@@ -0,0 +1,3 @@
+New: Added QGaussRadau quadrature up to and including 8 quadrature points. 
+<br>
+(Jan Philipp Thiele, 2023/09/15)
\ No newline at end of file
index e7ba3b2f851a5ca32a19a1a66b3178a25d27a02e..1ab6f439896c5d0c864e9038e9f421c797e66ddc 100644 (file)
@@ -47,6 +47,72 @@ public:
   QGauss(const unsigned int n);
 };
 
+/**
+ * The Gauss-Radau family of quadrature rules for numerical integration.
+ *
+ * This modification of the Gauss quadrature uses one of the two interval end
+ * points as well. Being exact for polynomials of degree $2n-2$, this
+ * formula is suboptimal by one degree.
+ *
+ * This formula is often used in the context of discontinuous Galerkin
+ * discretizations of ODEs and the temporal part of PDEs.
+ *
+ * The quadrature points are the left interval end point plus the $n-1$
+ * roots of the polynomial
+ * \f[
+ *   \frac{P_{n-1}(x)+P_n(x)}{1+x}
+ * \f]
+ * where $P_{n-1}$ and $P_n$ are Legendre polynomials.
+ * The quadrature weights are
+ * \f[
+ *   w_0=\frac{2}{n^2}\quad\text{and}
+ *   \quad w_i=\frac{1-x_i}{n^2(P_{n-1}(x_i))^2}\text{ for }i>0
+ * \f]
+ *
+ * For the right Gauss-Radau formula the quadrature points are
+ * $\tilde{x}_i=1-x_{n-i-1}$ and the weights are $\tilde{w}_i=w_{n-i-1}$,
+ * with $(x_i,w_i)$ as quadrature points
+ * and weights of the left Gauss-Radau formula.
+ *
+ * @see https://mathworld.wolfram.com/RadauQuadrature.html
+ */
+template <int dim>
+class QGaussRadau : public Quadrature<dim>
+{
+public:
+  /**
+   * EndPoint is used to specify which of the two endpoints of the unit interval
+   * is used also as quadrature point.
+   */
+  enum class EndPoint
+  {
+    /**
+     * Left end point.
+     */
+    left,
+    /**
+     * Right end point.
+     */
+    right
+  };
+  /**
+   * Generate a formula with <tt>n</tt> quadrature points (in each space
+   * direction).
+   * <tt>ep</tt> defines whether the left/lower/front endpoint(s) (default)
+   * or the right/upper/back endpoint(s) are part of the quadrature points.
+   */
+  QGaussRadau(const unsigned int n,
+              const EndPoint     end_point = QGaussRadau::EndPoint::left);
+
+  /**
+   * Move constructor.
+   */
+  QGaussRadau(QGaussRadau<dim> &&) noexcept = default;
+
+private:
+  const EndPoint end_point;
+};
+
 
 /**
  * The Gauss-Lobatto family of quadrature rules for numerical integration.
@@ -943,6 +1009,9 @@ public:
 template <>
 QGauss<1>::QGauss(const unsigned int n);
 template <>
+QGaussRadau<1>::QGaussRadau(const unsigned int             n,
+                            const QGaussRadau<1>::EndPoint end_poin);
+template <>
 QGaussLobatto<1>::QGaussLobatto(const unsigned int n);
 
 template <>
index 93bb3245365af65b15b6f9a1bf0d05e93edf3917..01242f3a61a42d395ce55c749ab925cf54e4e2e0 100644 (file)
@@ -144,6 +144,238 @@ namespace internal
 } // namespace internal
 
 
+namespace internal
+{
+  namespace QGaussRadau
+  {
+
+    // Implements lookup table after affine transformation to [0,1].
+    //
+    // Analytical values for [-1,1] and n < 4 listed on
+    // https://mathworld.wolfram.com/RadauQuadrature.html
+    // Values for n > 3 calculated with the Julia Package
+    // FastGaussQuadrature.jl
+    // https://github.com/JuliaApproximation/FastGaussQuadrature.jl
+    //
+    std::vector<double>
+    get_left_quadrature_points(const unsigned int n)
+    {
+      std::vector<double> q_points(n);
+      switch (n)
+        {
+          case 1:
+            q_points[0] = 0.;
+            break;
+          case 2:
+            q_points[0] = 0.;
+            q_points[1] = 2. / 3.;
+            break;
+          case 3:
+            q_points[0] = 0.;
+            q_points[1] = (6. - std::sqrt(6)) * 0.1;
+            q_points[2] = (6. + std::sqrt(6)) * 0.1;
+            break;
+
+          case 4:
+            q_points[0] = 0.000000000000000000;
+            q_points[1] = 0.212340538239152943;
+            q_points[2] = 0.590533135559265343;
+            q_points[3] = 0.911412040487296071;
+            break;
+          case 5:
+            q_points[0] = 0.000000000000000000;
+            q_points[1] = 0.139759864343780571;
+            q_points[2] = 0.416409567631083166;
+            q_points[3] = 0.723156986361876197;
+            q_points[4] = 0.942895803885482331;
+            break;
+          case 6:
+            q_points[0] = 0.000000000000000000;
+            q_points[1] = 0.098535085798826416;
+            q_points[2] = 0.304535726646363913;
+            q_points[3] = 0.562025189752613841;
+            q_points[4] = 0.801986582126391845;
+            q_points[5] = 0.960190142948531222;
+            break;
+          case 7:
+            q_points[0] = 0.000000000000000000;
+            q_points[1] = 0.073054328680258851;
+            q_points[2] = 0.230766137969945495;
+            q_points[3] = 0.441328481228449865;
+            q_points[4] = 0.663015309718845702;
+            q_points[5] = 0.851921400331515644;
+            q_points[6] = 0.970683572840215114;
+            break;
+          case 8:
+            q_points[0] = 0.000000000000000000;
+            q_points[1] = 0.056262560536922135;
+            q_points[2] = 0.180240691736892389;
+            q_points[3] = 0.352624717113169672;
+            q_points[4] = 0.547153626330555420;
+            q_points[5] = 0.734210177215410598;
+            q_points[6] = 0.885320946839095790;
+            q_points[7] = 0.977520613561287499;
+            break;
+          default:
+            Assert(false, ExcNotImplemented());
+            break;
+        }
+      return q_points;
+    }
+
+    std::vector<double>
+    get_quadrature_points(const unsigned int                       n,
+                          const ::dealii::QGaussRadau<1>::EndPoint end_point)
+    {
+      std::vector<double> left_points = get_left_quadrature_points(n);
+      switch (end_point)
+        {
+          case ::dealii::QGaussRadau<1>::EndPoint::left:
+            return left_points;
+          case ::dealii::QGaussRadau<1>::EndPoint::right:
+            {
+              std::vector<double> points(n);
+              for (unsigned int i = 0; i < n; ++i)
+                {
+                  points[n - i - 1] = 1. - left_points[i];
+                }
+              return points;
+            }
+          default:
+            Assert(
+              false,
+              ExcMessage(
+                "This constructor can only be called with either "
+                "QGaussRadau::left or QGaussRadau::right as second argument."));
+            return {};
+        }
+    }
+
+    // Implements lookup table after affine transformation to [0,1].
+    //
+    // Analytical values for [-1,1] and n < 4 listed on
+    // https://mathworld.wolfram.com/RadauQuadrature.html
+    // Values for n > 3 calculated with the Julia Package
+    // FastGaussQuadrature.jl
+    // https://github.com/JuliaApproximation/FastGaussQuadrature.jl
+    //
+    std::vector<double>
+    get_left_quadrature_weights(const unsigned int n)
+    {
+      std::vector<double> weights(n);
+      switch (n)
+        {
+          case 1:
+            weights[0] = 1.;
+            break;
+          case 2:
+            weights[0] = 0.25;
+            weights[1] = 0.75;
+            break;
+          case 3:
+            weights[0] = 1. / 9.;
+            weights[1] = (16. + std::sqrt(6)) / 36.;
+            weights[2] = (16. - std::sqrt(6)) / 36.;
+            break;
+          case 4:
+            weights[0] = 0.062500000000000000;
+            weights[1] = 0.328844319980059696;
+            weights[2] = 0.388193468843171852;
+            weights[3] = 0.220462211176768369;
+            break;
+          case 5:
+            weights[0] = 0.040000000000000001;
+            weights[1] = 0.223103901083570894;
+            weights[2] = 0.311826522975741427;
+            weights[3] = 0.281356015149462124;
+            weights[4] = 0.143713560791225797;
+            break;
+          case 6:
+            weights[0] = 0.027777777777777776;
+            weights[1] = 0.159820376610255471;
+            weights[2] = 0.242693594234484888;
+            weights[3] = 0.260463391594787597;
+            weights[4] = 0.208450667155953895;
+            weights[5] = 0.100794192626740456;
+            break;
+          case 7:
+            weights[0] = 0.020408163265306121;
+            weights[1] = 0.119613744612656100;
+            weights[2] = 0.190474936822115581;
+            weights[3] = 0.223554914507283209;
+            weights[4] = 0.212351889502977870;
+            weights[5] = 0.159102115733650767;
+            weights[6] = 0.074494235556010341;
+            break;
+          case 8:
+            weights[0] = 0.015625000000000000;
+            weights[1] = 0.092679077401489660;
+            weights[2] = 0.152065310323392683;
+            weights[3] = 0.188258772694559262;
+            weights[4] = 0.195786083726246729;
+            weights[5] = 0.173507397817250691;
+            weights[6] = 0.124823950664932445;
+            weights[7] = 0.057254407372128648;
+            break;
+
+          default:
+            Assert(false, dealii::StandardExceptions::ExcNotImplemented());
+            break;
+        }
+      return weights;
+    }
+
+    std::vector<double>
+    get_quadrature_weights(const unsigned int                       n,
+                           const ::dealii::QGaussRadau<1>::EndPoint end_point)
+    {
+      std::vector<double> left_weights = get_left_quadrature_weights(n);
+      switch (end_point)
+        {
+          case ::dealii::QGaussRadau<1>::EndPoint::left:
+            return left_weights;
+          case ::dealii::QGaussRadau<1>::EndPoint::right:
+            {
+              std::vector<double> weights(n);
+              for (unsigned int i = 0; i < n; ++i)
+                {
+                  weights[n - i - 1] = left_weights[i];
+                }
+              return weights;
+            }
+          default:
+            Assert(false,
+                   ExcMessage(
+                     "This constructor can only be called with either "
+                     "QGaussRadau::EndPoint::left or "
+                     "QGaussRadau::EndPoint::right as second argument."));
+            return {};
+        }
+    }
+  } // namespace QGaussRadau
+} // namespace internal
+
+#ifndef DOXYGEN
+template <>
+QGaussRadau<1>::QGaussRadau(const unsigned int n, const EndPoint end_point)
+  : Quadrature<1>(n)
+  , end_point(end_point)
+{
+  Assert(n > 0, ExcMessage("Need at least one point for quadrature rules."));
+  std::vector<double> p =
+    internal::QGaussRadau::get_quadrature_points(n, end_point);
+  std::vector<double> w =
+    internal::QGaussRadau::get_quadrature_weights(n, end_point);
+
+  for (unsigned int i = 0; i < this->size(); ++i)
+    {
+      this->quadrature_points[i] = dealii::Point<1>(p[i]);
+      this->weights[i]           = w[i];
+    }
+}
+#endif
+
+
 #ifndef DOXYGEN
 template <>
 QGaussLobatto<1>::QGaussLobatto(const unsigned int n)
@@ -820,6 +1052,15 @@ QGauss<dim>::QGauss(const unsigned int n)
 
 
 
+template <int dim>
+QGaussRadau<dim>::QGaussRadau(const unsigned int n, EndPoint end_point)
+  : Quadrature<dim>(
+      QGaussRadau<1>(n, static_cast<QGaussRadau<1>::EndPoint>(end_point)))
+  , end_point(end_point)
+{}
+
+
+
 template <int dim>
 QGaussLobatto<dim>::QGaussLobatto(const unsigned int n)
   : Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
@@ -2177,6 +2418,7 @@ QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
 // explicit specialization
 // note that 1d formulae are specialized by implementation above
 template class QGauss<2>;
+template class QGaussRadau<2>;
 template class QGaussLobatto<2>;
 template class QMidpoint<2>;
 template class QTrapezoid<2>;
@@ -2185,6 +2427,7 @@ template class QMilne<2>;
 template class QWeddle<2>;
 
 template class QGauss<3>;
+template class QGaussRadau<3>;
 template class QGaussLobatto<3>;
 template class QMidpoint<3>;
 template class QTrapezoid<3>;
index b190bb9b165c7a01c88b39344368c5f7791222c2..9aa68ee810f8bed6e387db0a72fe090b984e98d0 100644 (file)
@@ -46,6 +46,12 @@ fill_vector(std::vector<Quadrature<dim> *> &quadratures)
     {
       quadratures.push_back(new QGaussLobatto<dim>(i));
     }
+  for (unsigned int i = 1; i < 9; ++i)
+    {
+      quadratures.push_back(new QGaussRadau<dim>(i));
+      quadratures.push_back(
+        new QGaussRadau<dim>(i, QGaussRadau<dim>::EndPoint::right));
+    }
 }
 
 template <int dim>
index c93bc2d4cd8fa505f656258f515691733006d751..42383d7cf9b3824cec8ec7871e8fd228302811e6 100644 (file)
@@ -20,6 +20,22 @@ DEAL:1d::Quadrature no.17 is exact for polynomials of degree 5
 DEAL:1d::Quadrature no.18 is exact for polynomials of degree 7
 DEAL:1d::Quadrature no.19 is exact for polynomials of degree 9
 DEAL:1d::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:1d::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:1d::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:1d::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:1d::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:1d::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:1d::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:1d::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:1d::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:1d::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:1d::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:1d::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:1d::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:1d::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:1d::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:1d::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:1d::Quadrature no.36 is exact for polynomials of degree 14
 DEAL:2d::Quadrature no.0 is exact for polynomials of degree 1
 DEAL:2d::Quadrature no.1 is exact for polynomials of degree 1
 DEAL:2d::Quadrature no.2 is exact for polynomials of degree 3
@@ -41,6 +57,22 @@ DEAL:2d::Quadrature no.17 is exact for polynomials of degree 5
 DEAL:2d::Quadrature no.18 is exact for polynomials of degree 7
 DEAL:2d::Quadrature no.19 is exact for polynomials of degree 9
 DEAL:2d::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:2d::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:2d::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:2d::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:2d::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:2d::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:2d::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:2d::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:2d::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:2d::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:2d::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:2d::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:2d::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:2d::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:2d::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:2d::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:2d::Quadrature no.36 is exact for polynomials of degree 14
 DEAL:2d:faces::Quadrature no.0 is exact for polynomials of degree 1
 DEAL:2d:faces::Quadrature no.1 is exact for polynomials of degree 1
 DEAL:2d:faces::Quadrature no.2 is exact for polynomials of degree 3
@@ -62,6 +94,22 @@ DEAL:2d:faces::Quadrature no.17 is exact for polynomials of degree 5
 DEAL:2d:faces::Quadrature no.18 is exact for polynomials of degree 7
 DEAL:2d:faces::Quadrature no.19 is exact for polynomials of degree 9
 DEAL:2d:faces::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:2d:faces::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:2d:faces::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:2d:faces::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:2d:faces::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:2d:faces::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:2d:faces::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:2d:faces::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:2d:faces::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:2d:faces::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:2d:faces::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:2d:faces::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:2d:faces::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:2d:faces::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:2d:faces::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:2d:faces::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:2d:faces::Quadrature no.36 is exact for polynomials of degree 14
 DEAL:2d:subfaces::Quadrature no.0 is exact for polynomials of degree 1
 DEAL:2d:subfaces::Quadrature no.1 is exact for polynomials of degree 1
 DEAL:2d:subfaces::Quadrature no.2 is exact for polynomials of degree 3
@@ -83,6 +131,22 @@ DEAL:2d:subfaces::Quadrature no.17 is exact for polynomials of degree 5
 DEAL:2d:subfaces::Quadrature no.18 is exact for polynomials of degree 7
 DEAL:2d:subfaces::Quadrature no.19 is exact for polynomials of degree 9
 DEAL:2d:subfaces::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:2d:subfaces::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:2d:subfaces::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:2d:subfaces::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:2d:subfaces::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:2d:subfaces::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:2d:subfaces::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:2d:subfaces::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:2d:subfaces::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:2d:subfaces::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:2d:subfaces::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:2d:subfaces::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:2d:subfaces::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:2d:subfaces::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:2d:subfaces::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:2d:subfaces::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:2d:subfaces::Quadrature no.36 is exact for polynomials of degree 14
 DEAL:3d::Quadrature no.0 is exact for polynomials of degree 1
 DEAL:3d::Quadrature no.1 is exact for polynomials of degree 1
 DEAL:3d::Quadrature no.2 is exact for polynomials of degree 3
@@ -104,6 +168,22 @@ DEAL:3d::Quadrature no.17 is exact for polynomials of degree 5
 DEAL:3d::Quadrature no.18 is exact for polynomials of degree 7
 DEAL:3d::Quadrature no.19 is exact for polynomials of degree 9
 DEAL:3d::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:3d::Quadrature no.21 is exact for polynomials of degree 0
+DEAL:3d::Quadrature no.22 is exact for polynomials of degree 0
+DEAL:3d::Quadrature no.23 is exact for polynomials of degree 2
+DEAL:3d::Quadrature no.24 is exact for polynomials of degree 2
+DEAL:3d::Quadrature no.25 is exact for polynomials of degree 4
+DEAL:3d::Quadrature no.26 is exact for polynomials of degree 4
+DEAL:3d::Quadrature no.27 is exact for polynomials of degree 6
+DEAL:3d::Quadrature no.28 is exact for polynomials of degree 6
+DEAL:3d::Quadrature no.29 is exact for polynomials of degree 8
+DEAL:3d::Quadrature no.30 is exact for polynomials of degree 8
+DEAL:3d::Quadrature no.31 is exact for polynomials of degree 10
+DEAL:3d::Quadrature no.32 is exact for polynomials of degree 10
+DEAL:3d::Quadrature no.33 is exact for polynomials of degree 12
+DEAL:3d::Quadrature no.34 is exact for polynomials of degree 12
+DEAL:3d::Quadrature no.35 is exact for polynomials of degree 14
+DEAL:3d::Quadrature no.36 is exact for polynomials of degree 14
 DEAL:3d:faces::Quadrature no.0 is exact for polynomials of degree 1
 DEAL:3d:faces::Quadrature no.1 is exact for polynomials of degree 1
 DEAL:3d:faces::Quadrature no.2 is exact for polynomials of degree 3
@@ -125,6 +205,22 @@ DEAL:3d:faces::Quadrature no.17 is exact for polynomials of degree 5
 DEAL:3d:faces::Quadrature no.18 is exact for polynomials of degree 7
 DEAL:3d:faces::Quadrature no.19 is exact for polynomials of degree 9
 DEAL:3d:faces::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:3d:faces::Quadrature no.21 is exact for polynomials of degree 1
+DEAL:3d:faces::Quadrature no.22 is exact for polynomials of degree 1
+DEAL:3d:faces::Quadrature no.23 is exact for polynomials of degree 3
+DEAL:3d:faces::Quadrature no.24 is exact for polynomials of degree 3
+DEAL:3d:faces::Quadrature no.25 is exact for polynomials of degree 5
+DEAL:3d:faces::Quadrature no.26 is exact for polynomials of degree 5
+DEAL:3d:faces::Quadrature no.27 is exact for polynomials of degree 7
+DEAL:3d:faces::Quadrature no.28 is exact for polynomials of degree 7
+DEAL:3d:faces::Quadrature no.29 is exact for polynomials of degree 9
+DEAL:3d:faces::Quadrature no.30 is exact for polynomials of degree 9
+DEAL:3d:faces::Quadrature no.31 is exact for polynomials of degree 11
+DEAL:3d:faces::Quadrature no.32 is exact for polynomials of degree 11
+DEAL:3d:faces::Quadrature no.33 is exact for polynomials of degree 13
+DEAL:3d:faces::Quadrature no.34 is exact for polynomials of degree 13
+DEAL:3d:faces::Quadrature no.35 is exact for polynomials of degree 15
+DEAL:3d:faces::Quadrature no.36 is exact for polynomials of degree 15
 DEAL:3d:subfaces::Quadrature no.0 is exact for polynomials of degree 1
 DEAL:3d:subfaces::Quadrature no.1 is exact for polynomials of degree 1
 DEAL:3d:subfaces::Quadrature no.2 is exact for polynomials of degree 3
@@ -146,3 +242,19 @@ DEAL:3d:subfaces::Quadrature no.17 is exact for polynomials of degree 5
 DEAL:3d:subfaces::Quadrature no.18 is exact for polynomials of degree 7
 DEAL:3d:subfaces::Quadrature no.19 is exact for polynomials of degree 9
 DEAL:3d:subfaces::Quadrature no.20 is exact for polynomials of degree 11
+DEAL:3d:subfaces::Quadrature no.21 is exact for polynomials of degree 1
+DEAL:3d:subfaces::Quadrature no.22 is exact for polynomials of degree 1
+DEAL:3d:subfaces::Quadrature no.23 is exact for polynomials of degree 3
+DEAL:3d:subfaces::Quadrature no.24 is exact for polynomials of degree 3
+DEAL:3d:subfaces::Quadrature no.25 is exact for polynomials of degree 5
+DEAL:3d:subfaces::Quadrature no.26 is exact for polynomials of degree 5
+DEAL:3d:subfaces::Quadrature no.27 is exact for polynomials of degree 7
+DEAL:3d:subfaces::Quadrature no.28 is exact for polynomials of degree 7
+DEAL:3d:subfaces::Quadrature no.29 is exact for polynomials of degree 9
+DEAL:3d:subfaces::Quadrature no.30 is exact for polynomials of degree 9
+DEAL:3d:subfaces::Quadrature no.31 is exact for polynomials of degree 11
+DEAL:3d:subfaces::Quadrature no.32 is exact for polynomials of degree 11
+DEAL:3d:subfaces::Quadrature no.33 is exact for polynomials of degree 13
+DEAL:3d:subfaces::Quadrature no.34 is exact for polynomials of degree 13
+DEAL:3d:subfaces::Quadrature no.35 is exact for polynomials of degree 15
+DEAL:3d:subfaces::Quadrature no.36 is exact for polynomials of degree 15

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.