From 2fd90fe5524064a32cc1d7e1be34fb94e288be50 Mon Sep 17 00:00:00 2001
From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Sun, 25 Feb 2018 17:27:30 +0100
Subject: [PATCH] Move computation of QGaussRadauChebyshev support points into
 helper namespace

---
 include/deal.II/base/quadrature_lib.h |   8 --
 source/base/quadrature_lib.cc         | 103 +++++++++++++-------------
 2 files changed, 53 insertions(+), 58 deletions(-)

diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h
index 6576104b63..2ff37bb651 100644
--- a/include/deal.II/base/quadrature_lib.h
+++ b/include/deal.II/base/quadrature_lib.h
@@ -552,14 +552,6 @@ public:
 
 private:
   const EndPoint ep;
-  /// Computes the points of the quadrature formula.
-  static std::vector<double>
-  get_quadrature_points(const unsigned int n, EndPoint ep);
-
-  /// Computes the weights of the quadrature formula.
-  static std::vector<double>
-  get_quadrature_weights(const unsigned int n, EndPoint ep);
-
 };
 
 /**
diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc
index d766912061..1dec63ee94 100644
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -1228,65 +1228,68 @@ QGaussChebyshev<dim>::QGaussChebyshev (const unsigned int n)
 {}
 
 
-
-
-
-template <>
-std::vector<double>
-QGaussRadauChebyshev<1>::get_quadrature_points(const unsigned int n,
-                                               EndPoint ep)
+namespace internal
 {
+  namespace QGaussRadauChebyshev
+  {
+    // Computes the points of the quadrature formula.
+    std::vector<double>
+    get_quadrature_points(const unsigned int n,
+                          ::dealii::QGaussRadauChebyshev<1>::EndPoint ep)
+    {
 
-  std::vector<double> points(n);
-  // n point quadrature: index from 0 to n-1
-  for (unsigned short i=0; i<n; ++i)
-    // 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)
-      {
-      case QGaussRadauChebyshev::left:
-      {
-        points[i] = 1./2.*(1.-std::cos(numbers::PI*(1+2*double(i)/(2*double(n-1)+1.))));
-        break;
-      }
-
-      case QGaussRadauChebyshev::right:
-      {
-        points[i] = 1./2.*(1.-std::cos(numbers::PI*(2*double(n-1-i)/(2*double(n-1)+1.))));
-        break;
-      }
+      std::vector<double> points(n);
+      // n point quadrature: index from 0 to n-1
+      for (unsigned short i=0; i<n; ++i)
+        // 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)
+          {
+          case ::dealii::QGaussRadauChebyshev<1>::left:
+          {
+            points[i] = 1./2.*(1.-std::cos(numbers::PI*(1+2*double(i)/(2*double(n-1)+1.))));
+            break;
+          }
 
-      default:
-        Assert (false, ExcMessage ("This constructor can only be called with either "
-                                   "QGaussRadauChebyshev::left or QGaussRadauChebyshev::right as "
-                                   "second argument."));
-      }
+          case ::dealii::QGaussRadauChebyshev<1>::right:
+          {
+            points[i] = 1./2.*(1.-std::cos(numbers::PI*(2*double(n-1-i)/(2*double(n-1)+1.))));
+            break;
+          }
 
-  return points;
-}
+          default:
+            Assert (false, ExcMessage ("This constructor can only be called with either "
+                                       "QGaussRadauChebyshev::left or QGaussRadauChebyshev::right as "
+                                       "second argument."));
+          }
 
+      return points;
+    }
 
-template <>
-std::vector<double>
-QGaussRadauChebyshev<1>::get_quadrature_weights(const unsigned int n,
-                                                EndPoint ep)
-{
 
-  std::vector<double> weights(n);
 
-  for (unsigned short i=0; i<n; ++i)
+    // Computes the weights of the quadrature formula.
+    std::vector<double>
+    get_quadrature_weights(const unsigned int n,
+                           ::dealii::QGaussRadauChebyshev<1>::EndPoint ep)
     {
-      // same weights as on [-1,1]
-      weights[i] = 2.*numbers::PI/double(2*(n-1)+1.);
-      if (ep==left && i==0)
-        weights[i] /= 2.;
-      else if (ep==right && i==(n-1))
-        weights[i] /= 2.;
-    }
 
-  return weights;
+      std::vector<double> weights(n);
+
+      for (unsigned short i=0; i<n; ++i)
+        {
+          // same weights as on [-1,1]
+          weights[i] = 2.*numbers::PI/double(2*(n-1)+1.);
+          if (ep==::dealii::QGaussRadauChebyshev<1>::left && i==0)
+            weights[i] /= 2.;
+          else if (ep==::dealii::QGaussRadauChebyshev<1>::right && i==(n-1))
+            weights[i] /= 2.;
+        }
 
+      return weights;
+    }
+  }
 }
 
 
@@ -1299,8 +1302,8 @@ QGaussRadauChebyshev<1>::QGaussRadauChebyshev(const unsigned int n,
 {
 
   Assert(n>0,ExcMessage("Need at least one point for quadrature rules"));
-  std::vector<double> p=get_quadrature_points(n,ep);
-  std::vector<double> w=get_quadrature_weights(n,ep);
+  std::vector<double> p=internal::QGaussRadauChebyshev::get_quadrature_points(n,ep);
+  std::vector<double> w=internal::QGaussRadauChebyshev::get_quadrature_weights(n,ep);
 
   for (unsigned int i=0; i<this->size(); ++i)
     {
-- 
2.39.5