From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Thu, 16 Nov 2017 07:38:25 +0000 (+0100)
Subject: Add Sacado math header
X-Git-Tag: v9.0.0-rc1~764^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5473%2Fhead;p=dealii.git

Add Sacado math header

This adds some missing fundamental math functions for the supported
Sacado numbers.
---

diff --git a/doc/doxygen/headers/automatic_and_symbolic_differentiation.h b/doc/doxygen/headers/automatic_and_symbolic_differentiation.h
index 9b6195edb8..e54972b15d 100644
--- a/doc/doxygen/headers/automatic_and_symbolic_differentiation.h
+++ b/doc/doxygen/headers/automatic_and_symbolic_differentiation.h
@@ -73,6 +73,8 @@
  *   use Adol-C numbers.
  * - adolc_product_types.h: Defines some product and scalar types that allow the use of
  *   Adol-C numbers in conjunction with the Tensor and SymmetricTensor classes.
+ * - sacado_math.h: Extension of the sacado math operations that allow these numbers to be used
+ *   consistently throughout the library.
  * - sacado_number_types.h: Implementation of the internal classes that define how we
  *   use the supported Sacado numbers.
  * - sacado_product_types.h: Defines some product and scalar types that allow the use of
diff --git a/include/deal.II/differentiation/ad.h b/include/deal.II/differentiation/ad.h
index ccfbc517b6..74a9421407 100644
--- a/include/deal.II/differentiation/ad.h
+++ b/include/deal.II/differentiation/ad.h
@@ -27,6 +27,7 @@
 #include <deal.II/differentiation/ad/adolc_number_types.h>
 #include <deal.II/differentiation/ad/adolc_product_types.h>
 
+#include <deal.II/differentiation/ad/sacado_math.h>
 #include <deal.II/differentiation/ad/sacado_number_types.h>
 #include <deal.II/differentiation/ad/sacado_product_types.h>
 
diff --git a/include/deal.II/differentiation/ad/sacado_math.h b/include/deal.II/differentiation/ad/sacado_math.h
new file mode 100644
index 0000000000..3b2598713f
--- /dev/null
+++ b/include/deal.II/differentiation/ad/sacado_math.h
@@ -0,0 +1,96 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_differentiation_ad_sacado_math_h
+#define dealii_differentiation_ad_sacado_math_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#include <deal.II/base/numbers.h>
+#include <deal.II/differentiation/ad/ad_number_traits.h>
+#include <deal.II/differentiation/ad/sacado_number_types.h>
+
+
+#ifndef DOXYGEN
+
+/**
+ * Define some missing fundamental math functions
+ */
+namespace std
+{
+
+  /**
+   * Implementation of the error function for real-valued Sacado numbers.
+   */
+  template<typename ADNumberType, typename = typename std::enable_if<
+             dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value &&
+             dealii::Differentiation::AD::is_real_valued_ad_number<ADNumberType>::value
+             >::type>
+  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 = typename std::enable_if<
+             dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value
+             >::type>
+  inline ADNumberType erfc(const ADNumberType &x)
+  {
+    return 1.0 - std::erf(x);
+  }
+
+} // namespace std
+
+#endif // DOXYGEN
+
+#endif // DEAL_II_WITH_TRILINOS
+
+#endif