From 8dde01e6214947063c9c07f6edb39023e9154532 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 18 Jul 2005 19:01:46 +0000 Subject: [PATCH] Implement inversion of rank-4 tensors. git-svn-id: https://svn.dealii.org/trunk@11165 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/symmetric_tensor.h | 131 +++++++++++++++++++ 1 file changed, 131 insertions(+) diff --git a/deal.II/base/include/base/symmetric_tensor.h b/deal.II/base/include/base/symmetric_tensor.h index 09d2882e0b..bd2776ea3d 100644 --- a/deal.II/base/include/base/symmetric_tensor.h +++ b/deal.II/base/include/base/symmetric_tensor.h @@ -23,6 +23,7 @@ template class SymmetricTensor; template SymmetricTensor<2,dim> unit_symmetric_tensor (); template SymmetricTensor<4,dim> deviator_tensor (); template SymmetricTensor<4,dim> identity_tensor (); +template SymmetricTensor<4,dim> invert (const SymmetricTensor<4,dim> &); template double trace (const SymmetricTensor<2,dim2> &); template SymmetricTensor<2,dim> @@ -906,6 +907,9 @@ class SymmetricTensor template friend SymmetricTensor<4,dim2> identity_tensor (); + + template + friend SymmetricTensor<4,dim2> invert (const SymmetricTensor<4,dim2> &); }; @@ -2197,6 +2201,133 @@ identity_tensor () +/** + * Invert a symmetric rank-4 tensor. Since symmetric rank-4 tensors are + * mappings from and to symmetric rank-2 tensors, they can have an + * inverse. This function computes it, if it exists, for the case that the + * dimension equals 1. + * + * If a tensor is not invertible, then the result is unspecified, but will + * likely contain the results of a division by zero or a very small number at + * the very least. + * + * @relates SymmetricTensor + * @author Wolfgang Bangerth, 2005 + */ +template <> +inline +SymmetricTensor<4,1> +invert (const SymmetricTensor<4,1> &t) +{ + SymmetricTensor<4,1> tmp; + tmp.data[0][0] = 1./t.data[0][0]; + return tmp; +} + + + +/** + * Invert a symmetric rank-4 tensor. Since symmetric rank-4 tensors are + * mappings from and to symmetric rank-2 tensors, they can have an + * inverse. This function computes it, if it exists, for the case that the + * dimension equals 2. + * + * If a tensor is not invertible, then the result is unspecified, but will + * likely contain the results of a division by zero or a very small number at + * the very least. + * + * @relates SymmetricTensor + * @author Wolfgang Bangerth, 2005 + */ +template <> +inline +SymmetricTensor<4,2> +invert (const SymmetricTensor<4,2> &t) +{ + SymmetricTensor<4,2> tmp; + + // inverting this tensor is a little more + // complicated than necessary, since we + // store the data of 't' as a 3x3 matrix + // t.data, but the product between a rank-4 + // and a rank-2 tensor is really not the + // product between this matrix and the + // 3-vector of a rhs, but rather + // + // B.vec = t.data * mult * A.vec + // + // where mult is a 3x3 matrix with + // entries [[1,0,0],[0,1,0],[0,0,2]] to + // capture the fact that we need to add up + // both the c_ij12*a_12 and the c_ij21*a_21 + // terms + // + // in addition, in this scheme, the + // identity tensor has the matrix + // representation mult^-1. + // + // the inverse of 't' therefore has the + // matrix representation + // + // inv.data = mult^-1 * t.data^-1 * mult^-1 + // + // in order to compute it, let's first + // compute the inverse of t.data and put it + // into tmp.data; at the end of the + // function we then scale the last row and + // column of the inverse by 1/2, + // corresponding to the left and right + // multiplication with mult^-1 + const double t4 = t.data[0][0]*t.data[1][1], + t6 = t.data[0][0]*t.data[1][2], + t8 = t.data[0][1]*t.data[1][0], + t00 = t.data[0][2]*t.data[1][0], + t01 = t.data[0][1]*t.data[2][0], + t04 = t.data[0][2]*t.data[2][0], + t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-t8*t.data[2][2]+ + t00*t.data[2][1]+t01*t.data[1][2]-t04*t.data[1][1]); + tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07; + tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07; + tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07; + tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07; + tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07; + tmp.data[1][2] = -(t6-t00)*t07; + tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07; + tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07; + tmp.data[2][2] = (t4-t8)*t07; + + // scale last row and column as mentioned + // above + tmp.data[2][0] /= 2; + tmp.data[2][1] /= 2; + tmp.data[0][2] /= 2; + tmp.data[1][2] /= 2; + tmp.data[2][2] /= 4; + + return tmp; +} + + + +/** + * Invert a symmetric rank-4 tensor. Since symmetric rank-4 tensors are + * mappings from and to symmetric rank-2 tensors, they can have an + * inverse. This function computes it, if it exists, for the case that the + * dimension equals 3. + * + * If a tensor is not invertible, then the result is unspecified, but will + * likely contain the results of a division by zero or a very small number at + * the very least. + * + * @relates SymmetricTensor + * @author Wolfgang Bangerth, 2005 + */ +template <> +SymmetricTensor<4,2> +invert (const SymmetricTensor<4,2> &t); +// this function is implemented in the .cc file + + /** * Return the tensor of rank 4 that is the outer product of the two tensors * given as arguments, i.e. the result $T=t1 \otimes t2$ satisfies -- 2.39.5