From ac9fbba475368d72e0946ef4124609fc32d0e447 Mon Sep 17 00:00:00 2001
From: Marc Fehling <mafehling.git@gmail.com>
Date: Mon, 30 Nov 2020 18:37:04 -0700
Subject: [PATCH] Use boost implementation of legendre polynomials instead of
 gsl.

---
 cmake/checks/check_01_cxx_features.cmake | 22 +++++-
 include/deal.II/base/std_cxx17/cmath.h   | 86 +++++++++++++++++++++---
 source/fe/fe_series_legendre.cc          | 19 +-----
 3 files changed, 102 insertions(+), 25 deletions(-)

diff --git a/cmake/checks/check_01_cxx_features.cmake b/cmake/checks/check_01_cxx_features.cmake
index aca3421b23..556a22da3a 100644
--- a/cmake/checks/check_01_cxx_features.cmake
+++ b/cmake/checks/check_01_cxx_features.cmake
@@ -25,6 +25,7 @@
 #   DEAL_II_HAVE_FP_EXCEPTIONS
 #   DEAL_II_HAVE_COMPLEX_OPERATOR_OVERLOADS
 #   DEAL_II_HAVE_CXX17_BESSEL_FUNCTIONS
+#   DEAL_II_HAVE_CXX17_LEGENDRE_FUNCTIONS
 #   DEAL_II_FALLTHROUGH
 #   DEAL_II_DEPRECATED
 #   DEAL_II_CONSTEXPR
@@ -355,6 +356,7 @@ UNSET_IF_CHANGED(CHECK_CXX_FEATURES_FLAGS_SAVED
   DEAL_II_HAVE_CXX17_ATTRIBUTE_FALLTHROUGH
   DEAL_II_HAVE_ATTRIBUTE_FALLTHROUGH
   DEAL_II_HAVE_CXX17_BESSEL_FUNCTIONS
+  DEAL_II_HAVE_CXX17_LEGENDRE_FUNCTIONS
   DEAL_II_CXX14_CONSTEXPR_BUG_OK
   )
 
@@ -551,7 +553,7 @@ ENDIF()
 
 
 #
-# Check for c++17 bessel function support. Unfortunately libc++ version 10
+# Check for c++17 Bessel function support. Unfortunately libc++ version 10
 # does not have those.
 #
 
@@ -569,6 +571,24 @@ CHECK_CXX_SOURCE_COMPILES(
   )
 
 
+#
+# Check for c++17 Legendre function support.
+#
+
+CHECK_CXX_SOURCE_COMPILES(
+  "
+  #include <cmath>
+  using std::legendre;
+  using std::legendref;
+  using std::legendrel;
+  int main()
+  {
+  }
+  "
+  DEAL_II_HAVE_CXX17_LEGENDRE_FUNCTIONS
+  )
+
+
 #
 # Check for correct c++14 constexpr support.
 #
diff --git a/include/deal.II/base/std_cxx17/cmath.h b/include/deal.II/base/std_cxx17/cmath.h
index 6e402ac9ac..726a3a594e 100644
--- a/include/deal.II/base/std_cxx17/cmath.h
+++ b/include/deal.II/base/std_cxx17/cmath.h
@@ -17,38 +17,50 @@
 
 #include <deal.II/base/config.h>
 
-#ifdef DEAL_II_HAVE_CXX17_BESSEL_FUNCTIONS
+#if defined(DEAL_II_HAVE_CXX17_BESSEL_FUNCTIONS) || \
+  defined(DEAL_II_HAVE_CXX17_LEGENDRE_FUNCTIONS)
 #  include <cmath>
-#else
+#endif
+
+#ifndef DEAL_II_HAVE_CXX17_BESSEL_FUNCTIONS
 #  include <boost/math/special_functions/bessel.hpp>
 #endif
 
+#ifndef DEAL_II_HAVE_CXX17_LEGENDRE_FUNCTIONS
+#  include <deal.II/base/exceptions.h>
+
+#  include <boost/math/special_functions/legendre.hpp>
+
+#  include <limits>
+#endif
+
 
 DEAL_II_NAMESPACE_OPEN
+
 namespace std_cxx17
 {
 #ifndef DEAL_II_HAVE_CXX17_BESSEL_FUNCTIONS
 
   inline double
-  cyl_bessel_j(double x, double y)
+  cyl_bessel_j(double nu, double x)
   {
-    return boost::math::cyl_bessel_j(x, y);
+    return boost::math::cyl_bessel_j(nu, x);
   }
 
 
 
   inline float
-  cyl_bessel_jf(float x, float y)
+  cyl_bessel_jf(float nu, float x)
   {
-    return boost::math::cyl_bessel_j(x, y);
+    return boost::math::cyl_bessel_j(nu, x);
   }
 
 
 
   inline long double
-  cyl_bessel_jl(long double x, long double y)
+  cyl_bessel_jl(long double nu, long double x)
   {
-    return boost::math::cyl_bessel_j(x, y);
+    return boost::math::cyl_bessel_j(nu, x);
   }
 
 #else
@@ -56,7 +68,65 @@ namespace std_cxx17
   using std::cyl_bessel_jf;
   using std::cyl_bessel_jl;
 #endif
+
+#ifndef DEAL_II_HAVE_CXX17_LEGENDRE_FUNCTIONS
+
+  inline double
+  legendre(unsigned int l, double x)
+  {
+    Assert(static_cast<int>(l) >= 0,
+           ExcIndexRange(l, 0, std::numeric_limits<int>::max()));
+    return boost::math::legendre_p(static_cast<int>(l), x);
+  }
+
+
+
+  inline float
+  legendre(unsigned int l, float x)
+  {
+    Assert(static_cast<int>(l) >= 0,
+           ExcIndexRange(l, 0, std::numeric_limits<int>::max()));
+    return boost::math::legendre_p(static_cast<int>(l), x);
+  }
+
+
+
+  inline long double
+  legendre(unsigned int l, long double x)
+  {
+    Assert(static_cast<int>(l) >= 0,
+           ExcIndexRange(l, 0, std::numeric_limits<int>::max()));
+    return boost::math::legendre_p(static_cast<int>(l), x);
+  }
+
+
+
+  inline float
+  legendref(unsigned int l, float x)
+  {
+    Assert(static_cast<int>(l) >= 0,
+           ExcIndexRange(l, 0, std::numeric_limits<int>::max()));
+    return boost::math::legendre_p(static_cast<int>(l), x);
+  }
+
+
+
+  inline long double
+  legendrel(unsigned int l, long double x)
+  {
+    Assert(static_cast<int>(l) >= 0,
+           ExcIndexRange(l, 0, std::numeric_limits<int>::max()));
+    return boost::math::legendre_p(static_cast<int>(l), x);
+  }
+
+#else
+  using std::legendre;
+  using std::legendref;
+  using std::legendrel;
+#endif
 } // namespace std_cxx17
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif // dealii_cxx17_cmath_h
diff --git a/source/fe/fe_series_legendre.cc b/source/fe/fe_series_legendre.cc
index c2c63f0923..f74f217f3b 100644
--- a/source/fe/fe_series_legendre.cc
+++ b/source/fe/fe_series_legendre.cc
@@ -15,14 +15,12 @@
 
 
 
+#include <deal.II/base/std_cxx17/cmath.h>
 #include <deal.II/base/thread_management.h>
 
 #include <deal.II/fe/fe_series.h>
 
 #include <iostream>
-#ifdef DEAL_II_WITH_GSL
-#  include <gsl/gsl_sf_legendre.h>
-#endif
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -42,26 +40,15 @@ namespace
   double
   Lh(const Point<dim> &x_q, const TableIndices<dim> &indices)
   {
-#ifdef DEAL_II_WITH_GSL
     double res = 1.0;
     for (unsigned int d = 0; d < dim; d++)
       {
         const double x = 2.0 * (x_q[d] - 0.5);
         Assert((x_q[d] <= 1.0) && (x_q[d] >= 0.), ExcLegendre(d, x_q[d]));
-        const int ind = indices[d];
-        res *= std::sqrt(2.0) * gsl_sf_legendre_Pl(ind, x);
+        const unsigned int ind = indices[d];
+        res *= std::sqrt(2.0) * std_cxx17::legendre(ind, x);
       }
     return res;
-
-#else
-
-    (void)x_q;
-    (void)indices;
-    AssertThrow(false,
-                ExcMessage("deal.II has to be configured with GSL "
-                           "in order to use Legendre transformation."));
-    return 0;
-#endif
   }
 
 
-- 
2.39.5