From 5cb49a1adddcb6830a06eddebfe2948f616d9d97 Mon Sep 17 00:00:00 2001 From: Jan Philipp Thiele Date: Fri, 15 Sep 2023 11:30:17 +0200 Subject: [PATCH] Add QGaussRadau quadrature + tests --- doc/news/changes/minor/20230915Thiele | 3 + include/deal.II/base/quadrature_lib.h | 69 ++++++++ source/base/quadrature_lib.cc | 243 ++++++++++++++++++++++++++ tests/base/quadrature_test.cc | 6 + tests/base/quadrature_test.output | 112 ++++++++++++ 5 files changed, 433 insertions(+) create mode 100644 doc/news/changes/minor/20230915Thiele diff --git a/doc/news/changes/minor/20230915Thiele b/doc/news/changes/minor/20230915Thiele new file mode 100644 index 0000000000..6160f9be81 --- /dev/null +++ b/doc/news/changes/minor/20230915Thiele @@ -0,0 +1,3 @@ +New: Added QGaussRadau quadrature up to and including 8 quadrature points. +
+(Jan Philipp Thiele, 2023/09/15) \ No newline at end of file diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index e7ba3b2f85..1ab6f43989 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -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 +class QGaussRadau : public Quadrature +{ +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 n quadrature points (in each space + * direction). + * ep 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 &&) 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 <> diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 93bb324536..01242f3a61 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -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 + get_left_quadrature_points(const unsigned int n) + { + std::vector 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 + get_quadrature_points(const unsigned int n, + const ::dealii::QGaussRadau<1>::EndPoint end_point) + { + std::vector 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 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 + get_left_quadrature_weights(const unsigned int n) + { + std::vector 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 + get_quadrature_weights(const unsigned int n, + const ::dealii::QGaussRadau<1>::EndPoint end_point) + { + std::vector 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 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 p = + internal::QGaussRadau::get_quadrature_points(n, end_point); + std::vector 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::QGauss(const unsigned int n) +template +QGaussRadau::QGaussRadau(const unsigned int n, EndPoint end_point) + : Quadrature( + QGaussRadau<1>(n, static_cast::EndPoint>(end_point))) + , end_point(end_point) +{} + + + template QGaussLobatto::QGaussLobatto(const unsigned int n) : Quadrature(QGaussLobatto(n), QGaussLobatto<1>(n)) @@ -2177,6 +2418,7 @@ QGaussPyramid::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>; diff --git a/tests/base/quadrature_test.cc b/tests/base/quadrature_test.cc index b190bb9b16..9aa68ee810 100644 --- a/tests/base/quadrature_test.cc +++ b/tests/base/quadrature_test.cc @@ -46,6 +46,12 @@ fill_vector(std::vector *> &quadratures) { quadratures.push_back(new QGaussLobatto(i)); } + for (unsigned int i = 1; i < 9; ++i) + { + quadratures.push_back(new QGaussRadau(i)); + quadratures.push_back( + new QGaussRadau(i, QGaussRadau::EndPoint::right)); + } } template diff --git a/tests/base/quadrature_test.output b/tests/base/quadrature_test.output index c93bc2d4cd..42383d7cf9 100644 --- a/tests/base/quadrature_test.output +++ b/tests/base/quadrature_test.output @@ -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 -- 2.39.5