From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Thu, 12 Apr 2018 13:33:20 +0000 (+0200)
Subject: Fix Tensor::invert() for auto-differentiable numbers
X-Git-Tag: v9.0.0-rc1~180^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=35eb192fe928c176c36a6b9c113a6885acbe8672;p=dealii.git

Fix Tensor::invert() for auto-differentiable numbers
---

diff --git a/doc/news/changes/minor/20180412Jean-PaulPelteret-1 b/doc/news/changes/minor/20180412Jean-PaulPelteret-1
new file mode 100644
index 0000000000..0eadf9bd98
--- /dev/null
+++ b/doc/news/changes/minor/20180412Jean-PaulPelteret-1
@@ -0,0 +1,4 @@
+Fixed: The Tensor::invert() function would previously not work with some Sacado
+number types. This has now been fixed.
+<br>
+(Jean-Paul Pelteret, 2018/04/12)
diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index 6819b023f3..2aa612b6dc 100644
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -2065,7 +2065,7 @@ invert (const Tensor<2,1,Number> &t)
 {
   Number return_tensor [1][1];
 
-  return_tensor[0][0] = 1.0/t[0][0];
+  return_tensor[0][0] = internal::NumberType<Number>::value(1.0/t[0][0]);
 
   return Tensor<2,1,Number>(return_tensor);
 }
@@ -2080,7 +2080,7 @@ invert (const Tensor<2,2,Number> &t)
 
   // this is Maple output,
   // thus a bit unstructured
-  const Number inv_det_t = 1.0/(t[0][0]*t[1][1]-t[1][0]*t[0][1]);
+  const Number inv_det_t = internal::NumberType<Number>::value(1.0/(t[0][0]*t[1][1]-t[1][0]*t[0][1]));
   return_tensor[0][0] = t[1][1];
   return_tensor[0][1] = -t[0][1];
   return_tensor[1][0] = -t[1][0];
@@ -2098,23 +2098,24 @@ invert (const Tensor<2,3,Number> &t)
 {
   Tensor<2,3,Number> return_tensor;
 
-  const Number t4 = t[0][0]*t[1][1],
-               t6 = t[0][0]*t[1][2],
-               t8 = t[0][1]*t[1][0],
-               t00 = t[0][2]*t[1][0],
-               t01 = t[0][1]*t[2][0],
-               t04 = t[0][2]*t[2][0],
-               inv_det_t = 1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
-                                t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
-  return_tensor[0][0] = t[1][1]*t[2][2]-t[1][2]*t[2][1];
-  return_tensor[0][1] = t[0][2]*t[2][1]-t[0][1]*t[2][2];
-  return_tensor[0][2] = t[0][1]*t[1][2]-t[0][2]*t[1][1];
-  return_tensor[1][0] = t[1][2]*t[2][0]-t[1][0]*t[2][2];
-  return_tensor[1][1] = t[0][0]*t[2][2]-t04;
+  const Number t4 = internal::NumberType<Number>::value(t[0][0]*t[1][1]),
+               t6 = internal::NumberType<Number>::value(t[0][0]*t[1][2]),
+               t8 = internal::NumberType<Number>::value(t[0][1]*t[1][0]),
+               t00 = internal::NumberType<Number>::value(t[0][2]*t[1][0]),
+               t01 = internal::NumberType<Number>::value(t[0][1]*t[2][0]),
+               t04 = internal::NumberType<Number>::value(t[0][2]*t[2][0]),
+               inv_det_t = internal::NumberType<Number>::value(
+                             1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
+                                  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]));
+  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1]*t[2][2])-internal::NumberType<Number>::value(t[1][2]*t[2][1]);
+  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2]*t[2][1])-internal::NumberType<Number>::value(t[0][1]*t[2][2]);
+  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1]*t[1][2])-internal::NumberType<Number>::value(t[0][2]*t[1][1]);
+  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2]*t[2][0])-internal::NumberType<Number>::value(t[1][0]*t[2][2]);
+  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0]*t[2][2])-t04;
   return_tensor[1][2] = t00-t6;
-  return_tensor[2][0] = t[1][0]*t[2][1]-t[1][1]*t[2][0];
-  return_tensor[2][1] = t01-t[0][0]*t[2][1];
-  return_tensor[2][2] = t4-t8;
+  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0]*t[2][1])-internal::NumberType<Number>::value(t[1][1]*t[2][0]);
+  return_tensor[2][1] = t01-internal::NumberType<Number>::value(t[0][0]*t[2][1]);
+  return_tensor[2][2] = internal::NumberType<Number>::value(t4-t8);
   return_tensor *= inv_det_t;
 
   return return_tensor;