From 35eb192fe928c176c36a6b9c113a6885acbe8672 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Thu, 12 Apr 2018 15:33:20 +0200 Subject: [PATCH] Fix Tensor::invert() for auto-differentiable numbers --- .../changes/minor/20180412Jean-PaulPelteret-1 | 4 ++ include/deal.II/base/tensor.h | 37 ++++++++++--------- 2 files changed, 23 insertions(+), 18 deletions(-) create mode 100644 doc/news/changes/minor/20180412Jean-PaulPelteret-1 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. +
+(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::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::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::value(t[0][0]*t[1][1]), + t6 = internal::NumberType::value(t[0][0]*t[1][2]), + t8 = internal::NumberType::value(t[0][1]*t[1][0]), + t00 = internal::NumberType::value(t[0][2]*t[1][0]), + t01 = internal::NumberType::value(t[0][1]*t[2][0]), + t04 = internal::NumberType::value(t[0][2]*t[2][0]), + inv_det_t = internal::NumberType::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::value(t[1][1]*t[2][2])-internal::NumberType::value(t[1][2]*t[2][1]); + return_tensor[0][1] = internal::NumberType::value(t[0][2]*t[2][1])-internal::NumberType::value(t[0][1]*t[2][2]); + return_tensor[0][2] = internal::NumberType::value(t[0][1]*t[1][2])-internal::NumberType::value(t[0][2]*t[1][1]); + return_tensor[1][0] = internal::NumberType::value(t[1][2]*t[2][0])-internal::NumberType::value(t[1][0]*t[2][2]); + return_tensor[1][1] = internal::NumberType::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::value(t[1][0]*t[2][1])-internal::NumberType::value(t[1][1]*t[2][0]); + return_tensor[2][1] = t01-internal::NumberType::value(t[0][0]*t[2][1]); + return_tensor[2][2] = internal::NumberType::value(t4-t8); return_tensor *= inv_det_t; return return_tensor; -- 2.39.5