]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use a C++17 function instead of jn. 6996/head
authorDavid Wells <drwells@email.unc.edu>
Sat, 28 Jul 2018 17:00:35 +0000 (13:00 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sun, 29 Jul 2018 14:15:14 +0000 (10:15 -0400)
This is not guaranteed to be in the math library, but it is in C++17 (or in
boost).

cmake/checks/check_01_cxx_features.cmake
cmake/checks/check_02_system_features.cmake
include/deal.II/base/std_cxx17/cmath.h [new file with mode: 0644]
source/base/function_lib.cc
tests/base/functions_01.cc
tests/base/functions_01.output

index 8ef3144ec6d46b2c2d0b6e29336763c3fd79a140..32d61d3883c1ff9ca90546bb3a462391a43ef2bb 100644 (file)
@@ -176,6 +176,22 @@ IF(NOT DEFINED DEAL_II_WITH_CXX17 OR DEAL_II_WITH_CXX17)
       "
       DEAL_II_HAVE_CXX17_IF_CONSTEXPR)
 
+    #
+    # Not all compilers with C++17 support include the new special math functions:
+    #
+    CHECK_CXX_SOURCE_COMPILES(
+      "
+      #include <cmath>
+
+      int main()
+      {
+        std::cyl_bessel_j(1.0, 1.0);
+        std::cyl_bessel_jf(1.0f, 1.0f);
+        std::cyl_bessel_jl(1.0, 1.0);
+      }
+      "
+      DEAL_II_HAVE_CXX17_SPECIAL_MATH_FUNCTIONS)
+
     #
     # Some compilers treat lambdas as constexpr functions when compiling
     # with C++17 support even if they don't fulfill all the constexpr
@@ -211,6 +227,7 @@ IF(NOT DEFINED DEAL_II_WITH_CXX17 OR DEAL_II_WITH_CXX17)
 
   IF( DEAL_II_HAVE_CXX17_ATTRIBUTES AND
       DEAL_II_HAVE_CXX17_IF_CONSTEXPR AND
+      DEAL_II_HAVE_CXX17_SPECIAL_MATH_FUNCTIONS AND
       DEAL_II_NON_CONSTEXPR_LAMBDA)
     SET(DEAL_II_HAVE_CXX17 TRUE)
   ELSE()
index 6186752b43633fcce6f0fb10552b0697262ba004..4f0f832786e93dbc722295dd6b11bc4008c78cff 100644 (file)
@@ -18,7 +18,6 @@
 #
 #   DEAL_II_HAVE_GETHOSTNAME
 #   DEAL_II_HAVE_GETPID
-#   DEAL_II_HAVE_JN
 #   DEAL_II_HAVE_SYS_RESOURCE_H
 #   DEAL_II_HAVE_UNISTD_H
 #   DEAL_II_MSVC
@@ -40,23 +39,6 @@ CHECK_INCLUDE_FILE_CXX("unistd.h" DEAL_II_HAVE_UNISTD_H)
 CHECK_CXX_SYMBOL_EXISTS("gethostname" "unistd.h" DEAL_II_HAVE_GETHOSTNAME)
 CHECK_CXX_SYMBOL_EXISTS("getpid" "unistd.h" DEAL_II_HAVE_GETPID)
 
-#
-# Do we have the Bessel function jn?
-#
-FIND_SYSTEM_LIBRARY(m_LIBRARY NAMES m)
-MARK_AS_ADVANCED(m_LIBRARY)
-
-IF(NOT m_LIBRARY MATCHES "-NOTFOUND")
-  LIST(APPEND CMAKE_REQUIRED_LIBRARIES ${m_LIBRARY})
-  CHECK_CXX_SYMBOL_EXISTS("jn" "math.h" DEAL_II_HAVE_JN)
-  RESET_CMAKE_REQUIRED()
-  IF(DEAL_II_HAVE_JN)
-    LIST(APPEND DEAL_II_LIBRARIES ${m_LIBRARY})
-  ENDIF()
-ENDIF()
-
-
-
 ########################################################################
 #                                                                      #
 #                        Mac OSX specific setup:                       #
diff --git a/include/deal.II/base/std_cxx17/cmath.h b/include/deal.II/base/std_cxx17/cmath.h
new file mode 100644 (file)
index 0000000..1245a0b
--- /dev/null
@@ -0,0 +1,44 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+#ifndef dealii_cxx17_cmath_h
+#define dealii_cxx17_cmath_h
+
+#include <deal.II/base/config.h>
+
+#ifndef DEAL_II_WITH_CXX17
+#  include <boost/math/special_functions/bessel.hpp>
+#endif // DEAL_II_WITH_CXX17
+
+#include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+namespace std_cxx17
+{
+#ifndef DEAL_II_WITH_CXX17
+  double (&cyl_bessel_j)(double,
+                         double) = boost::math::cyl_bessel_j<double, double>;
+  float (&cyl_bessel_jf)(float,
+                         float)  = boost::math::cyl_bessel_j<float, float>;
+  long double (&cyl_bessel_jl)(long double, long double) =
+    boost::math::cyl_bessel_j<long double, long double>;
+#else
+  using std::cyl_bessel_j;
+  using std::cyl_bessel_jf;
+  using std::cyl_bessel_jl;
+#endif
+} // namespace std_cxx17
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // dealii_cxx17_cmath_h
index 02e6d621555bcffaff64fb41a4e0ad2cd14aa340..465cde311254944eacc8be7fa8d2b64bc77d8579 100644 (file)
@@ -17,6 +17,7 @@
 #include <deal.II/base/function_lib.h>
 #include <deal.II/base/numbers.h>
 #include <deal.II/base/point.h>
+#include <deal.II/base/std_cxx17/cmath.h>
 #include <deal.II/base/tensor.h>
 
 #include <deal.II/lac/vector.h>
@@ -2283,12 +2284,7 @@ namespace Functions
   {
     Assert(dim == 2, ExcNotImplemented());
     const double r = p.distance(center);
-#ifndef DEAL_II_HAVE_JN
-    Assert(false, ExcMessage("The Bessel function jn was not found by CMake."));
-    return r;
-#else
-    return jn(order, r * wave_number);
-#endif
+    return std_cxx17::cyl_bessel_j(order, r * wave_number);
   }
 
 
@@ -2300,17 +2296,11 @@ namespace Functions
   {
     Assert(dim == 2, ExcNotImplemented());
     AssertDimension(points.size(), values.size());
-#ifndef DEAL_II_HAVE_JN
-    (void)points;
-    (void)values;
-    Assert(false, ExcMessage("The Bessel function jn was not found by CMake."));
-#else
     for (unsigned int k = 0; k < points.size(); ++k)
       {
         const double r = points[k].distance(center);
-        values[k]      = jn(order, r * wave_number);
+        values[k]      = std_cxx17::cyl_bessel_j(order, r * wave_number);
       }
-#endif
   }
 
 
@@ -2319,23 +2309,19 @@ namespace Functions
   Bessel1<dim>::gradient(const Point<dim> &p, const unsigned int) const
   {
     Assert(dim == 2, ExcNotImplemented());
-#ifndef DEAL_II_HAVE_JN
-    (void)p;
-    Assert(false, ExcMessage("The Bessel function jn was not found by CMake."));
-    return Tensor<1, dim>();
-#else
     const double r  = p.distance(center);
     const double co = (r == 0.) ? 0. : (p(0) - center(0)) / r;
     const double si = (r == 0.) ? 0. : (p(1) - center(1)) / r;
 
-    const double dJn = (order == 0) ? (-jn(1, r * wave_number)) :
-                                      (.5 * (jn(order - 1, wave_number * r) -
-                                             jn(order + 1, wave_number * r)));
+    const double dJn =
+      (order == 0) ?
+        (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
+        (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
+               std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
     Tensor<1, dim> result;
     result[0] = wave_number * co * dJn;
     result[1] = wave_number * si * dJn;
     return result;
-#endif
   }
 
 
@@ -2348,9 +2334,6 @@ namespace Functions
   {
     Assert(dim == 2, ExcNotImplemented());
     AssertDimension(points.size(), gradients.size());
-#ifndef DEAL_II_HAVE_JN
-    Assert(false, ExcMessage("The Bessel function jn was not found by CMake."));
-#else
     for (unsigned int k = 0; k < points.size(); ++k)
       {
         const Point<dim> &p  = points[k];
@@ -2358,15 +2341,15 @@ namespace Functions
         const double      co = (r == 0.) ? 0. : (p(0) - center(0)) / r;
         const double      si = (r == 0.) ? 0. : (p(1) - center(1)) / r;
 
-        const double dJn = (order == 0) ?
-                             (-jn(1, r * wave_number)) :
-                             (.5 * (jn(order - 1, wave_number * r) -
-                                    jn(order + 1, wave_number * r)));
+        const double dJn =
+          (order == 0) ?
+            (-std_cxx17::cyl_bessel_j(1, r * wave_number)) :
+            (.5 * (std_cxx17::cyl_bessel_j(order - 1, wave_number * r) -
+                   std_cxx17::cyl_bessel_j(order + 1, wave_number * r)));
         Tensor<1, dim> &result = gradients[k];
         result[0]              = wave_number * co * dJn;
         result[1]              = wave_number * si * dJn;
       }
-#endif
   }
 
 
@@ -2837,6 +2820,7 @@ namespace Functions
   template class Monomial<1>;
   template class Monomial<2>;
   template class Monomial<3>;
+  template class Bessel1<1>;
   template class Bessel1<2>;
   template class Bessel1<3>;
   template class InterpolatedTensorProductGridData<1>;
index a3c269fdffb1b5866fae1bd42091537704d308e0..c0ca0e1dccc398b2eaa700765225e5e4fea32826 100644 (file)
 #include <deal.II/base/auto_derivative_function.h>
 #include <deal.II/base/data_out_base.h>
 #include <deal.II/base/flow_function.h>
+#include <deal.II/base/function_bessel.h>
 #include <deal.II/base/function_lib.h>
 #include <deal.II/base/job_identifier.h>
 #include <deal.II/base/quadrature_lib.h>
 
 #include <deal.II/lac/vector.h>
 
+#include <boost/core/demangle.hpp>
+
 #include <string>
+#include <typeinfo>
 #include <vector>
 
 #include "../tests.h"
 #include "functions.h"
 
-#define CHECK(F)                               \
-  {                                            \
-    deallog << #F << std::endl;                \
-    F f;                                       \
-    check_function_value_consistency(f, 5);    \
-    check_function_gradient_consistency(f, 5); \
-    check_gradient(f, 5);                      \
-  }
-
+template <typename FunctionType, typename... Args>
+void
+check(Args... args)
+{
+  FunctionType f(args...);
+  std::string  function_name = boost::core::demangle(typeid(f).name());
+  // remove the preceding dealii:: to match older output
+  function_name.erase(0, sizeof("dealii::") - 1 /*remember the NUL*/);
+  deallog << function_name << std::endl;
+  check_function_value_consistency(f, 5);
+  check_function_gradient_consistency(f, 5);
+  check_gradient(f, 5);
+}
 
-#define CHECKN(F, arg)                               \
-  {                                                  \
-    deallog << #F << '(' << arg << ')' << std::endl; \
-    F f(arg);                                        \
-    check_function_value_consistency(f, arg + 1);    \
-    check_function_gradient_consistency(f, arg + 1); \
-  }
 
 
 int
 main()
 {
-  std::string   logname = "output";
-  std::ofstream logfile(logname.c_str());
-  deallog.attach(logfile);
+  initlog();
 
-  CHECK(Functions::SquareFunction<1>);
-  CHECK(Functions::SquareFunction<2>);
-  CHECK(Functions::SquareFunction<3>);
-  CHECK(Functions::CosineFunction<1>);
-  CHECK(Functions::CosineFunction<2>);
-  CHECK(Functions::CosineFunction<3>);
-  CHECK(Functions::CosineGradFunction<1>);
-  CHECK(Functions::CosineGradFunction<2>);
-  CHECK(Functions::CosineGradFunction<3>);
-  CHECK(Functions::ExpFunction<1>);
-  CHECK(Functions::ExpFunction<2>);
-  CHECK(Functions::ExpFunction<3>);
+  check<Functions::SquareFunction<1>>();
+  check<Functions::SquareFunction<2>>();
+  check<Functions::SquareFunction<3>>();
+  check<Functions::CosineFunction<1>>();
+  check<Functions::CosineFunction<2>>();
+  check<Functions::CosineFunction<3>>();
+  check<Functions::CosineGradFunction<1>>();
+  check<Functions::CosineGradFunction<2>>();
+  check<Functions::CosineGradFunction<3>>();
+  check<Functions::ExpFunction<1>>();
+  check<Functions::ExpFunction<2>>();
+  check<Functions::ExpFunction<3>>();
+  // Bessel1 is only implemented in 2D right now
+  // check<Functions::Bessel1<1>>(3, 2.0, Point<1>(1.0));
+  check<Functions::Bessel1<2>>(3, 2.0, Point<2>(1.0, 2.0));
+  // check<Functions::Bessel1<3>>(3, 2.0, Point<2>(1.0, 2.0));
+  // check(Functions::Bessel1<3>);
 }
index 867ad3cb31b5210430fc96182474ddaf69b91761..ddca8cb87d86d7941e1238d93fd6ed56077b56bb 100644 (file)
@@ -71,3 +71,9 @@ DEAL::value list vs vector value list
 DEAL::gradient vs vector gradient list
 DEAL::gradient list vs vector gradient list
 DEAL::gradients vs difference quotients
+DEAL::Functions::Bessel1<2>
+DEAL::value vs vector value list
+DEAL::value list vs vector value list
+DEAL::gradient vs vector gradient list
+DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.