From a626e3706abbb13de110194f8de744b4f955b323 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 14 Apr 2000 12:01:59 +0000 Subject: [PATCH] Add a function that inverts tensors of rank 2. git-svn-id: https://svn.dealii.org/trunk@2720 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/tensor.h | 65 ++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/deal.II/base/include/base/tensor.h b/deal.II/base/include/base/tensor.h index 4fc6362ee1..83ddcba199 100644 --- a/deal.II/base/include/base/tensor.h +++ b/deal.II/base/include/base/tensor.h @@ -665,4 +665,69 @@ double determinant (const Tensor<2,3> &t) }; + +/** + * Compute and return the inverse of the given tensor. Since the + * compiler can perform the return value optimization, and since the + * size of the return object is known, it is acceptable to return the + * result by value, rather than by reference as a parameter. + * + * @author Wolfgang Bangerth, 2000 + */ +template +inline +Tensor<2,dim> +invert (const Tensor<2,dim> &t) +{ + Tensor<2,dim> return_tensor; + switch (dim) + { + case 1: + return_tensor[0][0] = 1.0/t[0][0]; + return return_tensor; + case 2: + // this is Maple output, + // thus a bit unstructured + { + const double t4 = 1.0/(t[0][0]*t[1][1]-t[0][1]*t[1][0]); + return_tensor[0][0] = t[1][1]*t4; + return_tensor[0][1] = -t[0][1]*t4; + return_tensor[1][0] = -t[1][0]*t4; + return_tensor[1][1] = t[0][0]*t4; + return return_tensor; + }; + + case 3: + { + const double 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], + t07 = 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])*t07; + return_tensor[0][1] = -(t[0][1]*t[2][2]-t[0][2]*t[2][1])*t07; + return_tensor[0][2] = -(-t[0][1]*t[1][2]+t[0][2]*t[1][1])*t07; + return_tensor[1][0] = -(t[1][0]*t[2][2]-t[1][2]*t[2][0])*t07; + return_tensor[1][1] = (t[0][0]*t[2][2]-t04)*t07; + return_tensor[1][2] = -(t6-t00)*t07; + return_tensor[2][0] = -(-t[1][0]*t[2][1]+t[1][1]*t[2][0])*t07; + return_tensor[2][1] = -(t[0][0]*t[2][1]-t01)*t07; + return_tensor[2][2] = (t4-t8)*t07; + return return_tensor; + }; + + // if desired, take over the + // inversion of a 4x4 tensor + // from the FullMatrix + + default: + AssertThrow (false, ExcNotImplemented()); + }; + return return_tensor; +}; + + #endif -- 2.39.5