From 6841cf1f7acab84a3624b4ff34c753a311844dd2 Mon Sep 17 00:00:00 2001 From: Jan Philipp Thiele Date: Mon, 18 Sep 2023 15:48:36 +0200 Subject: [PATCH] Unify QGaussRadauChebychev and QGaussRadau --- doc/news/changes/minor/20230916Thiele | 4 +++ include/deal.II/base/quadrature_lib.h | 19 +++++----- source/base/quadrature_lib.cc | 51 +++++++++++++++------------ tests/base/quadrature_chebyshev.cc | 3 +- 4 files changed, 44 insertions(+), 33 deletions(-) create mode 100644 doc/news/changes/minor/20230916Thiele diff --git a/doc/news/changes/minor/20230916Thiele b/doc/news/changes/minor/20230916Thiele new file mode 100644 index 0000000000..21010a5204 --- /dev/null +++ b/doc/news/changes/minor/20230916Thiele @@ -0,0 +1,4 @@ +Changed: The interface to QGaussRadauChebyshev +now matches that of the new quadrature QGaussRadau. +
+(Jan Philipp Thiele, 2023/09/16) diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 1ab6f43989..e0fbff7772 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -84,7 +84,7 @@ public: * EndPoint is used to specify which of the two endpoints of the unit interval * is used also as quadrature point. */ - enum class EndPoint + enum EndPoint { /** * Left end point. @@ -582,8 +582,9 @@ template class QGaussRadauChebyshev : public Quadrature { public: - /* EndPoint is used to specify which of the two endpoints of the unit interval - * is used also as quadrature point + /** + * EndPoint is used to specify which of the two endpoints of the unit interval + * is used also as quadrature point. */ enum EndPoint { @@ -597,17 +598,17 @@ public: right }; /// Generate a formula with n quadrature points - QGaussRadauChebyshev(const unsigned int n, - EndPoint ep = QGaussRadauChebyshev::left); + QGaussRadauChebyshev( + const unsigned int n, + const EndPoint end_point = QGaussRadauChebyshev::EndPoint::left); /** - * Move constructor. We cannot rely on the move constructor for `Quadrature`, - * since it does not know about the additional member `ep` of this class. + * Move constructor. */ QGaussRadauChebyshev(QGaussRadauChebyshev &&) noexcept = default; private: - const EndPoint ep; + const EndPoint end_point; }; /** @@ -1010,7 +1011,7 @@ template <> QGauss<1>::QGauss(const unsigned int n); template <> QGaussRadau<1>::QGaussRadau(const unsigned int n, - const QGaussRadau<1>::EndPoint end_poin); + const QGaussRadau<1>::EndPoint end_point); template <> QGaussLobatto<1>::QGaussLobatto(const unsigned int n); diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 01242f3a61..fd1207cff0 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1274,8 +1274,9 @@ namespace internal { // Computes the points of the quadrature formula. std::vector - get_quadrature_points(const unsigned int n, - ::dealii::QGaussRadauChebyshev<1>::EndPoint ep) + get_quadrature_points( + const unsigned int n, + const ::dealii::QGaussRadauChebyshev<1>::EndPoint end_point) { std::vector points(n); // n point quadrature: index from 0 to n-1 @@ -1283,9 +1284,9 @@ namespace internal // would be -cos(2i Pi/(2N+1)) // put + Pi so we start from the smallest point // then map from [-1,1] to [0,1] - switch (ep) + switch (end_point) { - case ::dealii::QGaussRadauChebyshev<1>::left: + case ::dealii::QGaussRadauChebyshev<1>::EndPoint::left: { points[i] = 1. / 2. * @@ -1295,7 +1296,7 @@ namespace internal break; } - case ::dealii::QGaussRadauChebyshev<1>::right: + case ::dealii::QGaussRadauChebyshev<1>::EndPoint::right: { points[i] = 1. / 2. * @@ -1309,8 +1310,8 @@ namespace internal false, ExcMessage( "This constructor can only be called with either " - "QGaussRadauChebyshev::left or QGaussRadauChebyshev::right as " - "second argument.")); + "QGaussRadauChebyshev::EndPoint::left or " + "QGaussRadauChebyshev::EndPoint:right as second argument.")); } return points; @@ -1320,8 +1321,9 @@ namespace internal // Computes the weights of the quadrature formula. std::vector - get_quadrature_weights(const unsigned int n, - ::dealii::QGaussRadauChebyshev<1>::EndPoint ep) + get_quadrature_weights( + const unsigned int n, + const ::dealii::QGaussRadauChebyshev<1>::EndPoint end_point) { std::vector weights(n); @@ -1329,9 +1331,11 @@ namespace internal { // same weights as on [-1,1] weights[i] = 2. * numbers::PI / double(2 * (n - 1) + 1.); - if (ep == ::dealii::QGaussRadauChebyshev<1>::left && i == 0) + if (end_point == ::dealii::QGaussRadauChebyshev<1>::EndPoint::left && + i == 0) weights[i] /= 2.; - else if (ep == ::dealii::QGaussRadauChebyshev<1>::right && + else if (end_point == + ::dealii::QGaussRadauChebyshev<1>::EndPoint::right && i == (n - 1)) weights[i] /= 2.; } @@ -1343,31 +1347,32 @@ namespace internal template <> -QGaussRadauChebyshev<1>::QGaussRadauChebyshev(const unsigned int n, EndPoint ep) +QGaussRadauChebyshev<1>::QGaussRadauChebyshev(const unsigned int n, + const EndPoint end_point) : Quadrature<1>(n) - , ep(ep) + , end_point(end_point) { - Assert(n > 0, ExcMessage("Need at least one point for quadrature rules")); - std::vector p = - internal::QGaussRadauChebyshev::get_quadrature_points(n, ep); - std::vector w = - internal::QGaussRadauChebyshev::get_quadrature_weights(n, ep); + Assert(n > 0, ExcMessage("Need at least one point for quadrature rules.")); + std::vector points = + internal::QGaussRadauChebyshev::get_quadrature_points(n, end_point); + std::vector weights = + internal::QGaussRadauChebyshev::get_quadrature_weights(n, end_point); for (unsigned int i = 0; i < this->size(); ++i) { - this->quadrature_points[i] = Point<1>(p[i]); - this->weights[i] = w[i]; + this->quadrature_points[i] = Point<1>(points[i]); + this->weights[i] = weights[i]; } } template QGaussRadauChebyshev::QGaussRadauChebyshev(const unsigned int n, - EndPoint ep) + EndPoint end_point) : Quadrature(QGaussRadauChebyshev<1>( n, - static_cast::EndPoint>(ep))) - , ep(ep) + static_cast::EndPoint>(end_point))) + , end_point(end_point) {} diff --git a/tests/base/quadrature_chebyshev.cc b/tests/base/quadrature_chebyshev.cc index 2e9e1c2690..ced9e24dd9 100644 --- a/tests/base/quadrature_chebyshev.cc +++ b/tests/base/quadrature_chebyshev.cc @@ -134,7 +134,8 @@ check_GRC_right(double *exact_monomials) { for (unsigned int n = 1; n < 18; ++n) { - QGaussRadauChebyshev<1> quadrature(n, QGaussRadauChebyshev<1>::right); + QGaussRadauChebyshev<1> quadrature( + n, QGaussRadauChebyshev<1>::EndPoint::right); const std::vector> &points = quadrature.get_points(); const std::vector &weights = quadrature.get_weights(); -- 2.39.5