From 47a70e19d08c683fd0f58a8a77406047b8e32d50 Mon Sep 17 00:00:00 2001
From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Fri, 27 Apr 2018 02:47:02 +0200
Subject: [PATCH] Adol-C advanced branching: Fix SymmetricTensor

---
 include/deal.II/base/symmetric_tensor.h       |  2 +-
 .../deal.II/base/symmetric_tensor.templates.h | 79 ++++++++++++++++++-
 2 files changed, 76 insertions(+), 5 deletions(-)

diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h
index d9611b895b..877bcc7cd5 100644
--- a/include/deal.II/base/symmetric_tensor.h
+++ b/include/deal.II/base/symmetric_tensor.h
@@ -1073,7 +1073,7 @@ inline
 SymmetricTensor<rank_,dim,Number> &
 SymmetricTensor<rank_,dim,Number>::operator = (const Number &d)
 {
-  Assert (d==internal::NumberType<Number>::value(0.0), ExcMessage ("Only assignment with zero is allowed"));
+  Assert (numbers::value_is_zero(d), ExcMessage ("Only assignment with zero is allowed"));
   (void) d;
 
   data = internal::NumberType<Number>::value(0.0);
diff --git a/include/deal.II/base/symmetric_tensor.templates.h b/include/deal.II/base/symmetric_tensor.templates.h
index 5d38eadbb2..80d5c7760d 100644
--- a/include/deal.II/base/symmetric_tensor.templates.h
+++ b/include/deal.II/base/symmetric_tensor.templates.h
@@ -42,7 +42,7 @@ std::array<Number,2>
 eigenvalues (const SymmetricTensor<2,2,Number> &T)
 {
   const Number upp_tri_sq = T[0][1]*T[0][1];
-  if (upp_tri_sq == Number(0.0))
+  if (upp_tri_sq == internal::NumberType<Number>::value(0.0))
     {
       // The tensor is diagonal
       std::array<Number,2> eig_vals =
@@ -59,7 +59,7 @@ eigenvalues (const SymmetricTensor<2,2,Number> &T)
       const Number tr_T = trace(T);
       const Number det_T = determinant(T);
       const Number descrim = tr_T*tr_T - 4.0*det_T;
-      Assert(descrim > Number(0.0), ExcMessage("The roots of the characteristic polynomial are complex valued."));
+      Assert(descrim > internal::NumberType<Number>::value(0.0), ExcMessage("The roots of the characteristic polynomial are complex valued."));
       const Number sqrt_desc = std::sqrt(descrim);
 
       const std::array<Number,2> eig_vals =
@@ -81,7 +81,7 @@ std::array<Number,3>
 eigenvalues (const SymmetricTensor<2,3,Number> &T)
 {
   const Number upp_tri_sq = T[0][1]*T[0][1] + T[0][2]*T[0][2] + T[1][2]*T[1][2];
-  if (upp_tri_sq == Number(0.0))
+  if (upp_tri_sq == internal::NumberType<Number>::value(0.0))
     {
       // The tensor is diagonal
       std::array<Number,3> eig_vals
@@ -711,7 +711,6 @@ namespace internal
     }
 
 
-
     template<typename Number>
     Tensor<2,1,Number>
     dediagonalize_tensor (const dealii::SymmetricTensor<2,1,Number> &T,
@@ -959,6 +958,78 @@ eigenvectors (const SymmetricTensor<2,dim,Number>         &T,
 
 
 
+#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
+
+// Specializations of the above functions for taped Adol-C numbers
+// when the advanced branching feature is activated.
+// We could copy-paste all of these functions and add the appropriate
+// conditional assignments (see the Adol-C manual, section 1.8).
+// However, some of the conditions are quite complicated (with possibly
+// no 1-1 correspondence for the operations along each branch) so this
+// needs some careful attention. For the sake of simplicity and until
+// this can be rigourously checked, at this point in time we simply
+// disable these functions.
+
+namespace internal
+{
+  namespace SymmetricTensorImplementation
+  {
+    template<>
+    struct Inverse<4,3,adouble>
+    {
+      static dealii::SymmetricTensor<4,3,adouble>
+      value (const dealii::SymmetricTensor<4,3,adouble> &/*t*/)
+      {
+        AssertThrow(false, ExcADOLCAdvancedBranching());
+        return dealii::SymmetricTensor<4,3,adouble>();
+      }
+    };
+  }
+}
+
+template<>
+std::array<adouble,1>
+eigenvalues (const SymmetricTensor<2,1,adouble> &/*T*/)
+{
+  AssertThrow(false, ExcADOLCAdvancedBranching());
+  return std::array<adouble,1> ();
+}
+
+
+
+template<>
+std::array<adouble,2>
+eigenvalues (const SymmetricTensor<2,2,adouble> &/*T*/)
+{
+  AssertThrow(false, ExcADOLCAdvancedBranching());
+  return std::array<adouble,2> ();
+}
+
+
+
+template<>
+std::array<adouble,3>
+eigenvalues (const SymmetricTensor<2,3,adouble> &/*T*/)
+{
+  AssertThrow(false, ExcADOLCAdvancedBranching());
+  return std::array<adouble,3> ();
+}
+
+
+
+template <int dim>
+std::array<std::pair<adouble, Tensor<1,dim,adouble> >,std::integral_constant<int, dim>::value>
+eigenvectors (const SymmetricTensor<2,dim,adouble>    &/*T*/,
+              const SymmetricTensorEigenvectorMethod   /*method*/)
+{
+  AssertThrow(false, ExcADOLCAdvancedBranching());
+  return std::array<std::pair<adouble, Tensor<1,dim,adouble> >,std::integral_constant<int, dim>::value> ();
+}
+
+#endif
+
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
-- 
2.39.5