From ac9fbba475368d72e0946ef4124609fc32d0e447 Mon Sep 17 00:00:00 2001 From: Marc Fehling 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 + 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 -#ifdef DEAL_II_HAVE_CXX17_BESSEL_FUNCTIONS +#if defined(DEAL_II_HAVE_CXX17_BESSEL_FUNCTIONS) || \ + defined(DEAL_II_HAVE_CXX17_LEGENDRE_FUNCTIONS) # include -#else +#endif + +#ifndef DEAL_II_HAVE_CXX17_BESSEL_FUNCTIONS # include #endif +#ifndef DEAL_II_HAVE_CXX17_LEGENDRE_FUNCTIONS +# include + +# include + +# include +#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(l) >= 0, + ExcIndexRange(l, 0, std::numeric_limits::max())); + return boost::math::legendre_p(static_cast(l), x); + } + + + + inline float + legendre(unsigned int l, float x) + { + Assert(static_cast(l) >= 0, + ExcIndexRange(l, 0, std::numeric_limits::max())); + return boost::math::legendre_p(static_cast(l), x); + } + + + + inline long double + legendre(unsigned int l, long double x) + { + Assert(static_cast(l) >= 0, + ExcIndexRange(l, 0, std::numeric_limits::max())); + return boost::math::legendre_p(static_cast(l), x); + } + + + + inline float + legendref(unsigned int l, float x) + { + Assert(static_cast(l) >= 0, + ExcIndexRange(l, 0, std::numeric_limits::max())); + return boost::math::legendre_p(static_cast(l), x); + } + + + + inline long double + legendrel(unsigned int l, long double x) + { + Assert(static_cast(l) >= 0, + ExcIndexRange(l, 0, std::numeric_limits::max())); + return boost::math::legendre_p(static_cast(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 #include #include #include -#ifdef DEAL_II_WITH_GSL -# include -#endif DEAL_II_NAMESPACE_OPEN @@ -42,26 +40,15 @@ namespace double Lh(const Point &x_q, const TableIndices &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