From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Fri, 27 Apr 2018 00:45:51 +0000 (+0200)
Subject: Add more helper functions to number.h
X-Git-Tag: v9.0.0-rc1~35^2~6
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=30555b7e27aa66ee09c9b0a1a08cd97429ccdef1;p=dealii.git

Add more helper functions to number.h
---

diff --git a/include/deal.II/base/numbers.h b/include/deal.II/base/numbers.h
index 542853c45a..a87894e737 100644
--- a/include/deal.II/base/numbers.h
+++ b/include/deal.II/base/numbers.h
@@ -192,6 +192,29 @@ namespace numbers
    */
   bool is_finite (const std::complex<long double> &x);
 
+  /**
+   * Return whether two numbers are equal to one another. For intricate data
+   * types (e.g. some automatically differentiable numbers), this function
+   * returns only whether the scalar values stored by the input values are
+   * equal.
+   *
+   * @note This function expects that @p value_2 is castable to the type
+   * of @p value_1.
+   */
+  template <typename Number1,typename Number2>
+  bool
+  values_are_equal (const Number1 &value_1, const Number2 &value_2);
+
+  /**
+   * Return whether or not a value is equal to zero. For intricate data
+   * types (e.g. some automatically differentiable numbers), this function
+   * returns only whether the scalar value stored by the input value is
+   * equal to zero.
+   */
+  template <typename Number>
+  bool
+  value_is_zero (const Number &value);
+
   /**
    * A structure that, together with its partial specializations
    * NumberTraits<std::complex<number> >, provides traits and member functions
@@ -542,6 +565,85 @@ namespace internal
   };
 }
 
+namespace numbers
+{
+
+#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
+
+  /**
+   * Return whether two numbers are equal to one another. For intricate data
+   * types (e.g. some automatically differentiable numbers), this function
+   * returns only whether the scalar values stored by the input values are
+   * equal.
+   *
+   * @note When Adol-C is compiled with the "advanced branching" feature, then
+   * this specialization is only intended for use in assertions and
+   * other code paths that do not affect the end result of a computation.
+   */
+  // Defined in differentiation/ad/adolc_number_types.cc
+  bool
+  values_are_equal (const adouble &value_1,
+                    const adouble &value_2);
+
+
+  /**
+   * Return whether two numbers are equal to one another. For intricate data
+   * types (e.g. some automatically differentiable numbers), this function
+   * returns only whether the scalar values stored by the input values are
+   * equal.
+   *
+   * @note When Adol-C is compiled with the "advanced branching" feature, then
+   * this specialization is only intended for use in assertions and
+   * other code paths that do not affect the end result of a computation.
+   */
+  template <typename Number>
+  bool
+  values_are_equal (const adouble &value_1,
+                    const Number  &value_2)
+  {
+    // Use the specialized definition for two Adol-C taped types
+    return values_are_equal(value_1, internal::NumberType<adouble>::value(value_2));
+  }
+
+
+  /**
+   * Return whether two numbers are equal to one another. For intricate data
+   * types (e.g. some automatically differentiable numbers), this function
+   * returns only whether the scalar values stored by the input values are
+   * equal.
+   *
+   * @note When Adol-C is compiled with the "advanced branching" feature, then
+   * this specialization is only intended for use in assertions and
+   * other code paths that do not affect the end result of a computation.
+   */
+  template <typename Number>
+  bool
+  values_are_equal (const Number  &value_1,
+                    const adouble &value_2)
+  {
+    // Use the above definition
+    return values_are_equal(value_2, value_1);
+  }
+
+#endif
+
+
+  template <typename Number1,typename Number2>
+  bool
+  values_are_equal (const Number1 &value_1, const Number2 &value_2)
+  {
+    return (value_1 == internal::NumberType<Number1>::value(value_2));
+  }
+
+
+  template <typename Number>
+  bool
+  value_is_zero (const Number &value)
+  {
+    return values_are_equal(value, 0.0);
+  }
+}
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
diff --git a/source/differentiation/ad/adolc_number_types.cc b/source/differentiation/ad/adolc_number_types.cc
index 073e421082..431ad64219 100644
--- a/source/differentiation/ad/adolc_number_types.cc
+++ b/source/differentiation/ad/adolc_number_types.cc
@@ -19,9 +19,55 @@
 #ifdef DEAL_II_WITH_ADOLC
 
 #include <deal.II/differentiation/ad/adolc_number_types.h>
+#include <deal.II/differentiation/ad/ad_number_traits.h>
+
+#include <utility>
 
 DEAL_II_NAMESPACE_OPEN
 
+
+#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
+
+namespace numbers
+{
+
+  namespace internal
+  {
+    namespace
+    {
+      // Apply some comparator and extract the boolean result of the operation,
+      // instead of the "adub" return type tpyically returned by Adol-C for
+      // such a comparison. This is implemented as a general function so that
+      // the list of implemented comparative operations can be easily extended.
+      bool
+      adouble_boolean_comparator (const adouble &value_1,
+                                  const adouble &value_2,
+                                  const std::function<adouble(const adouble &, const adouble &)> &comparator)
+      {
+        typedef typename Differentiation::AD::NumberTraits<double,Differentiation::AD::NumberTypes::adolc_taped>::ad_type ad_type;
+        static_assert(std::is_same<adouble,ad_type>::value, "The type of the AD number is not that which was expected.");
+        const ad_type result = comparator(value_1,value_2);
+        return !(Differentiation::AD::ADNumberTraits<ad_type>::get_scalar_value(result) == 0.0);
+      }
+    }
+  }
+
+  bool
+  values_are_equal (const adouble &value_1,
+                    const adouble &value_2)
+  {
+    return internal::adouble_boolean_comparator(
+             value_1, value_2,
+             [](const adouble &a,const adouble &b) -> adouble
+    {
+      return dealii::internal::NumberType<adouble>::value(a == b);
+    });
+  }
+}
+
+#endif
+
+
 /*---------------------- Explicit Instantiations ----------------------*/
 
 #include "adolc_number_types.inst"