"
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
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()
#
# 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
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: #
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
#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>
{
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);
}
{
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
}
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
}
{
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];
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
}
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>;
#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>);
}
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