Unify QGaussRadauChebychev and QGaussRadau

diff --git a/doc/news/changes/minor/20230916Thiele b/doc/news/changes/minor/20230916Thiele
new file mode 100644 (file)
index 0000000..21010a5
--- /dev/null
+++ b/doc/news/changes/minor/20230916Thiele
@@ -0,0 +1,4 @@
+Changed: The interface to QGaussRadauChebyshev<dim> 
+now matches that of the new quadrature QGaussRadau<dim>.
+<br>
+(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 1ab6f439896c5d0c864e9038e9f421c797e66ddc..e0fbff7772e97d2b12bc895f72f45b5e864c0932 100644 (file)
--- 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 <int dim>
 class QGaussRadauChebyshev : public Quadrature<dim>
 {
 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 <tt>n</tt> 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<dim> &&) 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 01242f3a61a42d395ce55c749ab925cf54e4e2e0..fd1207cff04404a7ff74f63fcfc9310551650eb4 100644 (file)
--- 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<double>
-    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<double> 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<double>
-    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<double> 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<double> p =
-    internal::QGaussRadauChebyshev::get_quadrature_points(n, ep);
-  std::vector<double> w =
-    internal::QGaussRadauChebyshev::get_quadrature_weights(n, ep);
+  Assert(n > 0, ExcMessage("Need at least one point for quadrature rules."));
+  std::vector<double> points =
+    internal::QGaussRadauChebyshev::get_quadrature_points(n, end_point);
+  std::vector<double> 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 <int dim>
 QGaussRadauChebyshev<dim>::QGaussRadauChebyshev(const unsigned int n,
-                                                EndPoint           ep)
+                                                EndPoint           end_point)
   : Quadrature<dim>(QGaussRadauChebyshev<1>(
       n,
-      static_cast<QGaussRadauChebyshev<1>::EndPoint>(ep)))
-  , ep(ep)
+      static_cast<QGaussRadauChebyshev<1>::EndPoint>(end_point)))
+  , end_point(end_point)
 {}
 
 

diff --git a/tests/base/quadrature_chebyshev.cc b/tests/base/quadrature_chebyshev.cc
index 2e9e1c26901653e40e2c98be9483144e7201b464..ced9e24dd9fff3cc32e30cc29d8306c96be5b7bf 100644 (file)
--- 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<Point<1>> &points  = quadrature.get_points();
       const std::vector<double>   &weights = quadrature.get_weights();