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

---
 include/deal.II/base/quadrature_lib.h |  54 -----
 source/base/quadrature_lib.cc         | 295 ++++++++++++++------------
 2 files changed, 161 insertions(+), 188 deletions(-)

diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h
index c330786217..fb373961e1 100644
--- a/include/deal.II/base/quadrature_lib.h
+++ b/include/deal.II/base/quadrature_lib.h
@@ -80,48 +80,6 @@ public:
    * direction).
    */
   QGaussLobatto(const unsigned int n);
-
-protected:
-  /**
-   * Compute Legendre-Gauss-Lobatto quadrature points in the interval $[-1,
-   * +1]$. They are equal to the roots of the corresponding Jacobi polynomial
-   * (specified by @p alpha, @p beta).  @p q is the number of points.
-   *
-   * @return Vector containing nodes.
-   */
-  std::vector<long double>
-  compute_quadrature_points (const unsigned int q,
-                             const int alpha,
-                             const int beta) const;
-
-  /**
-   * Compute Legendre-Gauss-Lobatto quadrature weights. The quadrature points
-   * and weights are related to Jacobi polynomial specified by @p alpha, @p
-   * beta. @p x denotes the quadrature points.
-   *
-   * @return Vector containing weights.
-   */
-  std::vector<long double>
-  compute_quadrature_weights (const std::vector<long double> &x,
-                              const int alpha,
-                              const int beta) const;
-
-  /**
-   * Evaluate a Jacobi polynomial $ P^{\alpha, \beta}_n(x) $ specified by the
-   * parameters @p alpha, @p beta, @p n. Note: The Jacobi polynomials are not
-   * orthonormal and defined on the interval $[-1, +1]$. @p x is the point of
-   * evaluation.
-   */
-  long double JacobiP(const long double x,
-                      const int alpha,
-                      const int beta,
-                      const unsigned int n) const;
-
-  /**
-   * Evaluate the Gamma function $ \Gamma(n) = (n-1)! $.
-   * @param n  point of evaluation (integer).
-   */
-  static long double gamma(const unsigned int n);
 };
 
 
@@ -878,18 +836,6 @@ public:
 
 template <> QGauss<1>::QGauss (const unsigned int n);
 template <> QGaussLobatto<1>::QGaussLobatto (const unsigned int n);
-template <>
-std::vector<long double> QGaussLobatto<1>::
-compute_quadrature_points(const unsigned int, const int, const int) const;
-template <>
-std::vector<long double> QGaussLobatto<1>::
-compute_quadrature_weights(const std::vector<long double> &, const int, const int) const;
-template <>
-long double QGaussLobatto<1>::
-JacobiP(const long double, const int, const int, const unsigned int) const;
-template <>
-long double
-QGaussLobatto<1>::gamma(const unsigned int n);
 
 template <> std::vector<double> QGaussLog<1>::get_quadrature_points(const unsigned int);
 template <> std::vector<double> QGaussLog<1>::get_quadrature_weights(const unsigned int);
diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc
index bf6c0c432e..45351f54b3 100644
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -161,169 +161,196 @@ QGauss<1>::QGauss (const unsigned int n)
     }
 }
 
-
-template <>
-QGaussLobatto<1>::QGaussLobatto (const unsigned int n)
-  :
-  Quadrature<1> (n)
+namespace internal
 {
-  Assert (n >= 2, ExcNotImplemented());
-
-  std::vector<long double> points  = compute_quadrature_points(n, 1, 1);
-  std::vector<long double> w       = compute_quadrature_weights(points, 0, 0);
-
-  // scale points to the interval
-  // [0.0, 1.0]:
-  for (unsigned int i=0; i<points.size(); ++i)
+  namespace QGaussLobatto
+  {
+    /*
+     * Evaluate a Jacobi polynomial $ P^{\alpha, \beta}_n(x) $ specified by the
+     * parameters @p alpha, @p beta, @p n. Note: The Jacobi polynomials are not
+     * orthonormal and defined on the interval $[-1, +1]$. @p x is the point of
+     * evaluation.
+     */
+    long double JacobiP(const long double x,
+                        const int alpha,
+                        const int beta,
+                        const unsigned int n)
     {
-      this->quadrature_points[i] = Point<1>(0.5 + 0.5*static_cast<double>(points[i]));
-      this->weights[i]           = 0.5*w[i];
-    }
-}
-
+      // the Jacobi polynomial is evaluated
+      // using a recursion formula.
+      std::vector<long double> p(n+1);
 
+      // initial values P_0(x), P_1(x):
+      p[0] = 1.0L;
+      if (n==0) return p[0];
+      p[1] = ((alpha+beta+2)*x + (alpha-beta))/2;
+      if (n==1) return p[1];
 
-template <>
-std::vector<long double> QGaussLobatto<1>::
-compute_quadrature_points(const unsigned int q,
-                          const int alpha,
-                          const int beta) const
-{
-  const unsigned int m = q-2; // no. of inner points
-  std::vector<long double> x(m);
-
-  // compute quadrature points with
-  // a Newton algorithm.
+      for (unsigned int i=1; i<=(n-1); ++i)
+        {
+          const int v  = 2*i + alpha + beta;
+          const int a1 = 2*(i+1)*(i + alpha + beta + 1)*v;
+          const int a2 = (v + 1)*(alpha*alpha - beta*beta);
+          const int a3 = v*(v + 1)*(v + 2);
+          const int a4 = 2*(i+alpha)*(i+beta)*(v + 2);
+
+          p[i+1] = static_cast<long double>( (a2 + a3*x)*p[i] - a4*p[i-1])/a1;
+        } // for
+      return p[n];
+    }
 
-  // Set tolerance. See class QGauss
-  // for detailed explanation.
-  const long double
-  long_double_eps = static_cast<long double>(std::numeric_limits<long double>::epsilon()),
-  double_eps      = static_cast<long double>(std::numeric_limits<double>::epsilon());
 
-  // check whether long double is
-  // more accurate than double, and
-  // set tolerances accordingly
-  volatile long double runtime_one = 1.0;
-  const long double tolerance
-    = (runtime_one + long_double_eps != runtime_one
-       ?
-       std::max (double_eps / 100, long_double_eps * 5)
-       :
-       double_eps * 5
-      );
 
-  // The following implementation
-  // follows closely the one given in
-  // the appendix of the book by
-  // Karniadakis and Sherwin:
-  // Spectral/hp element methods for
-  // computational fluid dynamics
-  // (Oxford University Press, 2005)
+    /**
+     * Evaluate the Gamma function $ \Gamma(n) = (n-1)! $.
+     * @param n  point of evaluation (integer).
+     */
+    long double gamma(const unsigned int n)
+    {
+      long double result = n - 1;
+      for (int i=n-2; i>1; --i)
+        result *= i;
+      return result;
+    }
 
-  // we take the zeros of the Chebyshev
-  // polynomial (alpha=beta=-0.5) as
-  // initial values:
-  for (unsigned int i=0; i<m; ++i)
-    x[i] = - std::cos( (long double) (2*i+1)/(2*m) * numbers::PI );
 
-  long double s, J_x, f, delta;
 
-  for (unsigned int k=0; k<m; ++k)
+    /**
+     * Compute Legendre-Gauss-Lobatto quadrature points in the interval $[-1,
+     * +1]$. They are equal to the roots of the corresponding Jacobi polynomial
+     * (specified by @p alpha, @p beta).  @p q is the number of points.
+     *
+     * @return Vector containing nodes.
+     */
+    std::vector<long double>
+    compute_quadrature_points(const unsigned int q,
+                              const int alpha,
+                              const int beta)
     {
-      long double r = x[k];
-      if (k>0)
-        r = (r + x[k-1])/2;
-
-      do
+      const unsigned int m = q-2; // no. of inner points
+      std::vector<long double> x(m);
+
+      // compute quadrature points with
+      // a Newton algorithm.
+
+      // Set tolerance. See class QGauss
+      // for detailed explanation.
+      const long double
+      long_double_eps = static_cast<long double>(std::numeric_limits<long double>::epsilon()),
+      double_eps      = static_cast<long double>(std::numeric_limits<double>::epsilon());
+
+      // check whether long double is
+      // more accurate than double, and
+      // set tolerances accordingly
+      volatile long double runtime_one = 1.0;
+      const long double tolerance
+        = (runtime_one + long_double_eps != runtime_one
+           ?
+           std::max (double_eps / 100, long_double_eps * 5)
+           :
+           double_eps * 5
+          );
+
+      // The following implementation
+      // follows closely the one given in
+      // the appendix of the book by
+      // Karniadakis and Sherwin:
+      // Spectral/hp element methods for
+      // computational fluid dynamics
+      // (Oxford University Press, 2005)
+
+      // we take the zeros of the Chebyshev
+      // polynomial (alpha=beta=-0.5) as
+      // initial values:
+      for (unsigned int i=0; i<m; ++i)
+        x[i] = - std::cos( (long double) (2*i+1)/(2*m) * numbers::PI );
+
+      long double s, J_x, f, delta;
+
+      for (unsigned int k=0; k<m; ++k)
         {
-          s = 0.;
-          for (unsigned int i=0; i<k; ++i)
-            s += 1./(r - x[i]);
-
-          J_x   =  0.5*(alpha + beta + m + 1)*JacobiP(r, alpha+1, beta+1, m-1);
-          f     = JacobiP(r, alpha, beta, m);
-          delta = f/(f*s- J_x);
-          r += delta;
-        }
-      while (std::fabs(delta) >= tolerance);
+          long double r = x[k];
+          if (k>0)
+            r = (r + x[k-1])/2;
 
-      x[k] = r;
-    } // for
+          do
+            {
+              s = 0.;
+              for (unsigned int i=0; i<k; ++i)
+                s += 1./(r - x[i]);
+
+              J_x   =  0.5*(alpha + beta + m + 1)*JacobiP(r, alpha+1, beta+1, m-1);
+              f     = JacobiP(r, alpha, beta, m);
+              delta = f/(f*s- J_x);
+              r += delta;
+            }
+          while (std::fabs(delta) >= tolerance);
 
-  // add boundary points:
-  x.insert(x.begin(), -1.L);
-  x.push_back(+1.L);
+          x[k] = r;
+        } // for
 
-  return x;
-}
+      // add boundary points:
+      x.insert(x.begin(), -1.L);
+      x.push_back(+1.L);
 
+      return x;
+    }
 
 
-template <>
-std::vector<long double> QGaussLobatto<1>::
-compute_quadrature_weights(const std::vector<long double> &x,
-                           const int alpha,
-                           const int beta) const
-{
-  const unsigned int q = x.size();
-  std::vector<long double> w(q);
-
-  const long double factor = std::pow(2., alpha+beta+1) *
-                             gamma(alpha+q) *
-                             gamma(beta+q) /
-                             ((q-1)*gamma(q)*gamma(alpha+beta+q+1));
-  for (unsigned int i=0; i<q; ++i)
+
+    /**
+     * Compute Legendre-Gauss-Lobatto quadrature weights. The quadrature points
+     * and weights are related to Jacobi polynomial specified by @p alpha, @p
+     * beta. @p x denotes the quadrature points.
+     *
+     * @return Vector containing weights.
+     */
+    std::vector<long double>
+    compute_quadrature_weights(const std::vector<long double> &x,
+                               const int alpha,
+                               const int beta)
     {
-      const long double s = JacobiP(x[i], alpha, beta, q-1);
-      w[i] = factor/(s*s);
-    }
-  w[0]   *= (beta + 1);
-  w[q-1] *= (alpha + 1);
+      const unsigned int q = x.size();
+      std::vector<long double> w(q);
+
+      const long double factor = std::pow(2., alpha+beta+1) *
+                                 gamma(alpha+q) *
+                                 gamma(beta+q) /
+                                 ((q-1)*gamma(q)*gamma(alpha+beta+q+1));
+      for (unsigned int i=0; i<q; ++i)
+        {
+          const long double s = JacobiP(x[i], alpha, beta, q-1);
+          w[i] = factor/(s*s);
+        }
+      w[0]   *= (beta + 1);
+      w[q-1] *= (alpha + 1);
 
-  return w;
+      return w;
+    }
+  }
 }
 
 
 
 template <>
-long double QGaussLobatto<1>::JacobiP(const long double x,
-                                      const int alpha,
-                                      const int beta,
-                                      const unsigned int n) const
+QGaussLobatto<1>::QGaussLobatto (const unsigned int n)
+  :
+  Quadrature<1> (n)
 {
-  // the Jacobi polynomial is evaluated
-  // using a recursion formula.
-  std::vector<long double> p(n+1);
+  Assert (n >= 2, ExcNotImplemented());
 
-  // initial values P_0(x), P_1(x):
-  p[0] = 1.0L;
-  if (n==0) return p[0];
-  p[1] = ((alpha+beta+2)*x + (alpha-beta))/2;
-  if (n==1) return p[1];
+  std::vector<long double> points
+    = internal::QGaussLobatto::compute_quadrature_points(n, 1, 1);
+  std::vector<long double> w
+    = internal::QGaussLobatto::compute_quadrature_weights(points, 0, 0);
 
-  for (unsigned int i=1; i<=(n-1); ++i)
+  // scale points to the interval
+  // [0.0, 1.0]:
+  for (unsigned int i=0; i<points.size(); ++i)
     {
-      const int v  = 2*i + alpha + beta;
-      const int a1 = 2*(i+1)*(i + alpha + beta + 1)*v;
-      const int a2 = (v + 1)*(alpha*alpha - beta*beta);
-      const int a3 = v*(v + 1)*(v + 2);
-      const int a4 = 2*(i+alpha)*(i+beta)*(v + 2);
-
-      p[i+1] = static_cast<long double>( (a2 + a3*x)*p[i] - a4*p[i-1])/a1;
-    } // for
-  return p[n];
-}
-
-
-
-template <>
-long double QGaussLobatto<1>::gamma(const unsigned int n)
-{
-  long double result = n - 1;
-  for (int i=n-2; i>1; --i)
-    result *= i;
-  return result;
+      this->quadrature_points[i] = Point<1>(0.5 + 0.5*static_cast<double>(points[i]));
+      this->weights[i]           = 0.5*w[i];
+    }
 }
 
 
-- 
2.39.5