From: Jean-Paul Pelteret Date: Wed, 16 Aug 2017 09:06:46 +0000 (-0600) Subject: Generalised dim 3 rank 4 symmetric tensor invert() function X-Git-Tag: v9.0.0-rc1~1196^2~8 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=98fe90dc44fb2ec5b059d5011a4ee46c631f2ce2;p=dealii.git Generalised dim 3 rank 4 symmetric tensor invert() function --- diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index 510fc2c91f..575f9d48f6 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -52,6 +52,20 @@ determinant (const SymmetricTensor<2,dim,Number> &); namespace internal { + /** + * A namespace for functions and classes that are internal to how the + * SymmetricTensor class (and its associate functions) works. + */ + namespace SymmetricTensor + { + /** + * Compute the inverse of a symmetric tensor of a + * generic @p rank, @p dim and @p Number type. + */ + template + struct Inverse; + } + /** * A namespace for classes that are internal to how the SymmetricTensor * class works. @@ -837,11 +851,13 @@ private: template friend SymmetricTensor<4,dim2,Number2> identity_tensor (); - template - friend SymmetricTensor<2,dim2,Number2> invert (const SymmetricTensor<2,dim2,Number2> &); - template - friend SymmetricTensor<4,dim2,Number2> invert (const SymmetricTensor<4,dim2,Number2> &); + /** + * Make a few helper classes friends as well. + */ + friend struct internal::SymmetricTensor::Inverse<2,dim,Number>; + + friend struct internal::SymmetricTensor::Inverse<4,dim,Number>; }; @@ -1070,6 +1086,271 @@ namespace internal return t; } + + + template + struct Inverse<2,1,Number> + { + static inline dealii::SymmetricTensor<2,1,Number> + value (const dealii::SymmetricTensor<2,1,Number> &t) + { + dealii::SymmetricTensor<2,1,Number> tmp; + + tmp[0][0] = 1.0/t[0][0]; + + return tmp; + } + }; + + + template + struct Inverse<2,2,Number> + { + static inline dealii::SymmetricTensor<2,2,Number> + value (const dealii::SymmetricTensor<2,2,Number> &t) + { + dealii::SymmetricTensor<2,2,Number> tmp; + + // Sympy result: ([ + // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)], + // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ]) + const TableIndices<2> idx_00 (0,0); + const TableIndices<2> idx_01 (0,1); + const TableIndices<2> idx_11 (1,1); + const Number inv_det_t + = 1.0/(t[idx_00]*t[idx_11] + - t[idx_01]*t[idx_01]); + tmp[idx_00] = t[idx_11]; + tmp[idx_01] = -t[idx_01]; + tmp[idx_11] = t[idx_00]; + tmp *= inv_det_t; + + return tmp; + } + }; + + + template + struct Inverse<2,3,Number> + { + static dealii::SymmetricTensor<2,3,Number> + value (const dealii::SymmetricTensor<2,3,Number> &t) + { + dealii::SymmetricTensor<2,3,Number> tmp; + + // Sympy result: ([ + // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), + // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), + // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)], + // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), + // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), + // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)], + // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), + // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11), + // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ]) + const TableIndices<2> idx_00 (0,0); + const TableIndices<2> idx_01 (0,1); + const TableIndices<2> idx_02 (0,2); + const TableIndices<2> idx_11 (1,1); + const TableIndices<2> idx_12 (1,2); + const TableIndices<2> idx_22 (2,2); + const Number inv_det_t + = 1.0/(t[idx_00]*t[idx_11]*t[idx_22] + - t[idx_00]*t[idx_12]*t[idx_12] + - t[idx_01]*t[idx_01]*t[idx_22] + + 2.0*t[idx_01]*t[idx_02]*t[idx_12] + - t[idx_02]*t[idx_02]*t[idx_11]); + tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12]; + tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12]; + tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11]; + tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02]; + tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02]; + tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01]; + tmp *= inv_det_t; + + return tmp; + } + }; + + + template + struct Inverse<4,1,Number> + { + static inline dealii::SymmetricTensor<4,1,Number> + value (const dealii::SymmetricTensor<4,1,Number> &t) + { + dealii::SymmetricTensor<4,1,Number> tmp; + tmp.data[0][0] = 1.0/t.data[0][0]; + return tmp; + } + }; + + + template + struct Inverse<4,2,Number> + { + static inline dealii::SymmetricTensor<4,2,Number> + value (const dealii::SymmetricTensor<4,2,Number> &t) + { + dealii::SymmetricTensor<4,2,Number> 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 Number 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; + } + }; + + + template + struct Inverse<4,3,Number> + { + static dealii::SymmetricTensor<4,3,Number> + value (const dealii::SymmetricTensor<4,3,Number> &t) + { + dealii::SymmetricTensor<4,3,Number> tmp = t; + + // This function follows the exact same scheme as the 2d case, except + // that hardcoding the inverse of a 6x6 matrix is pretty wasteful. + // Instead, we use the Gauss-Jordan algorithm implemented for + // FullMatrix. For historical reasons the following code is copied from + // there, with the tangential benefit that we do not need to copy the + // tensor entries to and from the FullMatrix. + const unsigned int N = 6; + + // First get an estimate of the size of the elements of this matrix, + // for later checks whether the pivot element is large enough, or + // whether we have to fear that the matrix is not regular. + Number diagonal_sum = internal::NumberType::value(0.0); + for (unsigned int i=0; i max) + { + max = std::fabs(tmp.data[i][j]); + r = i; + } + + // Check whether the pivot is too small + Assert(max > 1.e-16*typical_diagonal_element, + ExcMessage("This tensor seems to be noninvertible")); + + // Row interchange + if (r>j) + { + for (unsigned int k=0; k &T); namespace internal { - /** - * A namespace for functions and classes that are internal to how the - * SymmetricTensor class (and its associate functions) works. - */ namespace SymmetricTensor { /** @@ -3063,109 +3340,16 @@ identity_tensor () template inline SymmetricTensor<2,dim,Number> -invert (const SymmetricTensor<2,dim,Number> &) +invert (const SymmetricTensor<2,dim,Number> &t) { - // if desired, take over the - // inversion of a 4x4 tensor - // from the FullMatrix - AssertThrow (false, ExcNotImplemented()); - - return SymmetricTensor<2,dim,Number>(); + return internal::SymmetricTensor::Inverse<2,dim,Number>::value(t); } -#ifndef DOXYGEN - -template -inline -SymmetricTensor<2,1,Number> -invert (const SymmetricTensor<2,1,Number> &t) -{ - SymmetricTensor<2,1,Number> tmp; - - tmp[0][0] = 1.0/t[0][0]; - - return tmp; -} - - - -template -inline -SymmetricTensor<2,2,Number> -invert (const SymmetricTensor<2,2,Number> &t) -{ - SymmetricTensor<2,2,Number> tmp; - - // Sympy result: ([ - // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)], - // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ]) - const TableIndices<2> idx_00 (0,0); - const TableIndices<2> idx_01 (0,1); - const TableIndices<2> idx_11 (1,1); - const Number inv_det_t - = 1.0/(t[idx_00]*t[idx_11] - - t[idx_01]*t[idx_01]); - tmp[idx_00] = t[idx_11]; - tmp[idx_01] = -t[idx_01]; - tmp[idx_11] = t[idx_00]; - tmp *= inv_det_t; - - return tmp; -} - - - -template -inline -SymmetricTensor<2,3,Number> -invert (const SymmetricTensor<2,3,Number> &t) -{ - SymmetricTensor<2,3,Number> tmp; - - // Sympy result: ([ - // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), - // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), - // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)], - // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), - // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), - // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)], - // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11), - // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11), - // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ]) - const TableIndices<2> idx_00 (0,0); - const TableIndices<2> idx_01 (0,1); - const TableIndices<2> idx_02 (0,2); - const TableIndices<2> idx_11 (1,1); - const TableIndices<2> idx_12 (1,2); - const TableIndices<2> idx_22 (2,2); - const Number inv_det_t - = 1.0/(t[idx_00]*t[idx_11]*t[idx_22] - - t[idx_00]*t[idx_12]*t[idx_12] - - t[idx_01]*t[idx_01]*t[idx_22] - + 2.0*t[idx_01]*t[idx_02]*t[idx_12] - - t[idx_02]*t[idx_02]*t[idx_11]); - tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12]; - tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12]; - tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11]; - tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02]; - tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02]; - tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01]; - tmp *= inv_det_t; - - return tmp; -} - -#endif /* DOXYGEN */ - - - /** * 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 either 1 or 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 @@ -3179,103 +3363,11 @@ inline SymmetricTensor<4,dim,Number> invert (const SymmetricTensor<4,dim,Number> &t) { - SymmetricTensor<4,dim,Number> tmp; - switch (dim) - { - case 1: - tmp.data[0][0] = 1./t.data[0][0]; - break; - case 2: - - // 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 Number 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; - } - break; - default: - Assert (false, ExcNotImplemented()); - } - return tmp; + return internal::SymmetricTensor::Inverse<4,dim,Number>::value(t); } -/** - * 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,3,double> -invert (const SymmetricTensor<4,3,double> &t); -// this function is implemented in the .cc file for double data types - - - /** * 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 T phi = diff --git a/source/base/symmetric_tensor.cc b/source/base/symmetric_tensor.cc index cb167583aa..de0f31c6d1 100644 --- a/source/base/symmetric_tensor.cc +++ b/source/base/symmetric_tensor.cc @@ -31,116 +31,6 @@ DEAL_II_ENABLE_EXTRA_DIAGNOSTICS DEAL_II_NAMESPACE_OPEN -template <> -SymmetricTensor<4,3,double> -invert<3,double> (const SymmetricTensor<4,3,double> &t) -{ - SymmetricTensor<4,3,double> tmp = t; - - // this function follows the exact same - // scheme as the 2d case, except that - // hardcoding the inverse of a 6x6 matrix - // is pretty wasteful. instead, we use the - // Gauss-Jordan algorithm implemented for - // FullMatrix; the following code is copied - // from there because using the FullMatrix - // class would introduce circular - // references between libbase and liblac - const unsigned int N = 6; - - // first get an estimate of the - // size of the elements of this - // matrix, for later checks whether - // the pivot element is large - // enough, or whether we have to - // fear that the matrix is not - // regular - double diagonal_sum = 0; - for (unsigned int i=0; i max) - { - max = std::fabs(tmp.data[i][j]); - r = i; - } - // check whether the pivot is - // too small - Assert(max > 1.e-16*typical_diagonal_element, - ExcMessage("This tensor seems to be noninvertible")); - - // row interchange - if (r>j) - { - for (unsigned int k=0; k const unsigned int SymmetricTensor::dimension; @@ -149,6 +39,7 @@ template const unsigned int SymmetricTensor::n_independent_components; +// explicit instantiations #include "symmetric_tensor.inst" diff --git a/source/base/symmetric_tensor.inst.in b/source/base/symmetric_tensor.inst.in index f428a5516d..72679648cc 100644 --- a/source/base/symmetric_tensor.inst.in +++ b/source/base/symmetric_tensor.inst.in @@ -79,6 +79,19 @@ for (deal_II_dimension : DIMENSIONS; number : COMPLEX_SCALARS) for (number : REAL_SCALARS) { + namespace internal + \{ + namespace SymmetricTensor + \{ + template + struct Inverse<4,3,number>; + \} + \} + + template + SymmetricTensor<4,3,number> + invert (const SymmetricTensor<4,3,number> &t); + template std::array eigenvalues (const SymmetricTensor<2,1,number> &);