From a626e3706abbb13de110194f8de744b4f955b323 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
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 <int dim>
+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