# include <adolc/adouble.h> // Taped double
#endif
-// Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
-// but header indirectly references numbers.h. We therefore simply
-// import the whole Sacado header at this point to get the math
-// functions imported into the standard namespace.
-#ifdef DEAL_II_TRILINOS_WITH_SACADO
-DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
-# include <Sacado.hpp>
-DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
-#endif
namespace std
{
typename NumberTraits<number>::real_type
NumberTraits<number>::abs(const number &x)
{
- return std::abs(x);
+ // Make things work with AD types
+ using std::abs;
+ return abs(x);
}
typename NumberTraits<std::complex<number>>::real_type
NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
{
- return std::abs(x);
+ // Make things work with AD types
+ using std::abs;
+ return abs(x);
}
// whether we have to fear that the matrix is not regular.
Number diagonal_sum = internal::NumberType<Number>::value(0.0);
for (unsigned int i = 0; i < N; ++i)
- diagonal_sum += std::fabs(tmp.data[i][i]);
+ diagonal_sum += numbers::NumberTraits<Number>::abs(tmp.data[i][i]);
const Number typical_diagonal_element =
diagonal_sum / static_cast<double>(N);
(void)typical_diagonal_element;
{
// Pivot search: search that part of the line on and right of the
// diagonal for the largest element.
- Number max = std::fabs(tmp.data[j][j]);
- unsigned int r = j;
+ Number max = numbers::NumberTraits<Number>::abs(tmp.data[j][j]);
+ unsigned int r = j;
for (unsigned int i = j + 1; i < N; ++i)
- if (std::fabs(tmp.data[i][j]) > max)
+ if (numbers::NumberTraits<Number>::abs(tmp.data[i][j]) > max)
{
- max = std::fabs(tmp.data[i][j]);
+ max = numbers::NumberTraits<Number>::abs(tmp.data[i][j]);
r = i;
}
compute_norm(const typename SymmetricTensorAccessors::
StorageType<2, dim, Number>::base_tensor_type &data)
{
+ // Make things work with AD types
+ using std::sqrt;
switch (dim)
{
case 1:
return numbers::NumberTraits<Number>::abs(data[0]);
case 2:
- return std::sqrt(
- numbers::NumberTraits<Number>::abs_square(data[0]) +
- numbers::NumberTraits<Number>::abs_square(data[1]) +
- 2. * numbers::NumberTraits<Number>::abs_square(data[2]));
+ return sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
+ numbers::NumberTraits<Number>::abs_square(data[1]) +
+ 2. * numbers::NumberTraits<Number>::abs_square(data[2]));
case 3:
- return std::sqrt(
- numbers::NumberTraits<Number>::abs_square(data[0]) +
- numbers::NumberTraits<Number>::abs_square(data[1]) +
- numbers::NumberTraits<Number>::abs_square(data[2]) +
- 2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
- 2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
- 2. * numbers::NumberTraits<Number>::abs_square(data[5]));
+ return sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
+ numbers::NumberTraits<Number>::abs_square(data[1]) +
+ numbers::NumberTraits<Number>::abs_square(data[2]) +
+ 2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
+ 2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
+ 2. * numbers::NumberTraits<Number>::abs_square(data[5]));
default:
{
return_value +=
2. * numbers::NumberTraits<Number>::abs_square(data[d]);
- return std::sqrt(return_value);
+ return sqrt(return_value);
}
}
}
compute_norm(const typename SymmetricTensorAccessors::
StorageType<4, dim, Number>::base_tensor_type &data)
{
+ // Make things work with AD types
+ using std::sqrt;
switch (dim)
{
case 1:
return_value +=
4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
- return std::sqrt(return_value);
+ return sqrt(return_value);
}
}
}
inline typename numbers::NumberTraits<Number>::real_type
Tensor<rank_, dim, Number>::norm() const
{
- return std::sqrt(norm_square());
+ // Make things work with AD types
+ using std::sqrt;
+ return sqrt(norm_square());
}
{
Number sum = internal::NumberType<Number>::value(0.0);
for (unsigned int i = 0; i < dim; ++i)
- sum += std::fabs(t[i][j]);
+ sum += numbers::NumberTraits<Number>::abs(t[i][j]);
if (sum > max)
max = sum;
{
Number sum = internal::NumberType<Number>::value(0.0);
for (unsigned int j = 0; j < dim; ++j)
- sum += std::fabs(t[i][j]);
+ sum += numbers::NumberTraits<Number>::abs(t[i][j]);
if (sum > max)
max = sum;
# ifndef DOXYGEN
+DEAL_II_NAMESPACE_OPEN
/**
- * Define some missing fundamental math functions
+ * Implementation of the error function for real-valued Sacado numbers.
*/
-namespace std
+template <
+ typename ADNumberType,
+ typename = std::enable_if_t<
+ dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value &&
+ dealii::Differentiation::AD::is_real_valued_ad_number<ADNumberType>::value>>
+inline ADNumberType
+erf(ADNumberType x)
{
- /**
- * Implementation of the error function for real-valued Sacado numbers.
- */
- template <
- typename ADNumberType,
- typename = std::enable_if_t<
- dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value &&
- dealii::Differentiation::AD::is_real_valued_ad_number<
- ADNumberType>::value>>
- inline ADNumberType
- erf(ADNumberType x)
- {
- // Reference:
- // Handbook of Mathematical Functions: with Formulas, Graphs, and
- // Mathematical Tables Abramowitz, M. and Stegun, I. Dover Books on
- // Mathematics. 1972.
- //
- // Current implementation source:
- // https://www.johndcook.com/blog/cpp_erf/
- // https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
- //
- // Note: This implementation has a reported maximum round-off error
- // of 1.5e-7.
-
- // Constants
- const double a1 = 0.254829592;
- const double a2 = -0.284496736;
- const double a3 = 1.421413741;
- const double a4 = -1.453152027;
- const double a5 = 1.061405429;
- const double p = 0.3275911;
-
- // Save the sign of x
- const bool neg_val =
- (x < dealii::internal::NumberType<ADNumberType>::value(0.0) ? true :
- false);
- x = std::fabs(x);
-
- // Abramowitz1972a equation 7.1.26
- const ADNumberType t = 1.0 / (1.0 + p * x);
- const ADNumberType y =
- 1.0 -
- (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x);
-
- if (!neg_val)
- return y;
- else
- return -y;
- }
-
-
- /**
- * Implementation of the complementary error function for Sacado numbers.
- */
- template <
- typename ADNumberType,
- typename = std::enable_if_t<
- dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value>>
- inline ADNumberType
- erfc(const ADNumberType &x)
- {
- return 1.0 - std::erf(x);
- }
-
-} // namespace std
+ // Reference:
+ // Handbook of Mathematical Functions: with Formulas, Graphs, and
+ // Mathematical Tables Abramowitz, M. and Stegun, I. Dover Books on
+ // Mathematics. 1972.
+ //
+ // Current implementation source:
+ // https://www.johndcook.com/blog/cpp_erf/
+ // https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
+ //
+ // Note: This implementation has a reported maximum round-off error
+ // of 1.5e-7.
+
+ // Constants
+ const double a1 = 0.254829592;
+ const double a2 = -0.284496736;
+ const double a3 = 1.421413741;
+ const double a4 = -1.453152027;
+ const double a5 = 1.061405429;
+ const double p = 0.3275911;
+
+ // Save the sign of x
+ const bool neg_val =
+ (x < dealii::internal::NumberType<ADNumberType>::value(0.0) ? true : false);
+ x = std::fabs(x);
+
+ // Abramowitz1972a equation 7.1.26
+ const ADNumberType t = 1.0 / (1.0 + p * x);
+ const ADNumberType y =
+ 1.0 -
+ (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x);
+
+ if (!neg_val)
+ return y;
+ else
+ return -y;
+}
+
+
+/**
+ * Implementation of the complementary error function for Sacado numbers.
+ */
+template <typename ADNumberType,
+ typename = std::enable_if_t<
+ dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value>>
+inline ADNumberType
+erfc(const ADNumberType &x)
+{
+ return 1.0 - std::erf(x);
+}
+
+DEAL_II_NAMESPACE_CLOSE
# endif // DOXYGEN
projects / dealii.git / commitdiff