From: Nils Much Date: Tue, 15 Dec 2020 07:38:41 +0000 (+0100) Subject: Move and rename tensor basis transformations X-Git-Tag: v9.3.0-rc1~634^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6176f9aefb8a24513717349c9ad75e2d19ae7364;p=dealii.git Move and rename tensor basis transformations Add short documentation to basis transformation methods --- diff --git a/doc/news/changes/minor/20210112Much b/doc/news/changes/minor/20210112Much new file mode 100644 index 0000000000..51d7808d76 --- /dev/null +++ b/doc/news/changes/minor/20210112Much @@ -0,0 +1,5 @@ +New: The old tensor basis transformation functions internal::Physics::transformation_contraction() +have been moved out of the internal namespace and renamed to +Physics::Transformations::basis_transformation() and have documentation now. +
+(Nils Much, 2021/01/12) diff --git a/include/deal.II/physics/transformations.h b/include/deal.II/physics/transformations.h index d0a365c6f1..483f22ed3b 100644 --- a/include/deal.II/physics/transformations.h +++ b/include/deal.II/physics/transformations.h @@ -178,7 +178,7 @@ namespace Physics /** * Return the result of the push forward transformation on a rank-4 - * contravariant tensor, i.e. (in index notation) + * contravariant tensor, i.e. (in index notation): * @f[ * \left[ \chi\left(\bullet\right)^{\sharp} \right]_{ijkl} * \dealcoloneq F_{iI} F_{jJ} @@ -197,7 +197,7 @@ namespace Physics /** * Return the result of the push forward transformation on a rank-4 - * contravariant symmetric tensor, i.e. (in index notation) + * contravariant symmetric tensor, i.e. (in index notation): * @f[ * \left[ \chi\left(\bullet\right)^{\sharp} \right]_{ijkl} * \dealcoloneq F_{iI} F_{jJ} @@ -280,7 +280,7 @@ namespace Physics /** * Return the result of the pull back transformation on a rank-4 - * contravariant tensor, i.e. (in index notation) + * contravariant tensor, i.e. (in index notation): * @f[ * \left[ \chi^{-1}\left(\bullet\right)^{\sharp} \right]_{IJKL} * \dealcoloneq F^{-1}_{Ii} F^{-1}_{Jj} @@ -299,7 +299,7 @@ namespace Physics /** * Return the result of the pull back transformation on a rank-4 - * contravariant symmetric tensor, i.e. (in index notation) + * contravariant symmetric tensor, i.e. (in index notation): * @f[ * \left[ \chi^{-1}\left(\bullet\right)^{\sharp} \right]_{IJKL} * \dealcoloneq F^{-1}_{Ii} F^{-1}_{Jj} @@ -403,7 +403,7 @@ namespace Physics /** * Return the result of the push forward transformation on a rank-4 - * covariant tensor, i.e. (in index notation) + * covariant tensor, i.e. (in index notation): * @f[ * \left[ \chi\left(\bullet\right)^{\flat} \right]_{ijkl} * \dealcoloneq F^{-T}_{iI} F^{-T}_{jJ} @@ -422,7 +422,7 @@ namespace Physics /** * Return the result of the push forward transformation on a rank-4 - * covariant symmetric tensor, i.e. (in index notation) + * covariant symmetric tensor, i.e. (in index notation): * @f[ * \left[ \chi\left(\bullet\right)^{\flat} \right]_{ijkl} * \dealcoloneq F^{-T}_{iI} F^{-T}_{jJ} @@ -505,7 +505,7 @@ namespace Physics /** * Return the result of the pull back transformation on a rank-4 - * contravariant tensor, i.e. (in index notation) + * contravariant tensor, i.e. (in index notation): * @f[ * \left[ \chi^{-1}\left(\bullet\right)^{\flat} \right]_{IJKL} * \dealcoloneq F^{T}_{Ii} F^{T}_{Jj} @@ -524,7 +524,7 @@ namespace Physics /** * Return the result of the pull back transformation on a rank-4 - * contravariant symmetric tensor, i.e. (in index notation) + * contravariant symmetric tensor, i.e. (in index notation): * @f[ * \left[ \chi^{-1}\left(\bullet\right)^{\flat} \right]_{IJKL} * \dealcoloneq F^{T}_{Ii} F^{T}_{Jj} @@ -619,7 +619,7 @@ namespace Physics /** * Return the result of the push forward transformation on a rank-4 - * contravariant tensor, i.e. (in index notation) + * contravariant tensor, i.e. (in index notation): * @f[ * \textrm{det} \mathbf{F}^{-1} \; \left[ * \chi\left(\bullet\right)^{\sharp} \right]_{ijkl} @@ -640,7 +640,7 @@ namespace Physics /** * Return the result of the push forward transformation on a rank-4 - * contravariant symmetric tensor, i.e. (in index notation) + * contravariant symmetric tensor, i.e. (in index notation): * @f[ * \textrm{det} \mathbf{F}^{-1} \; \left[ * \chi\left(\bullet\right)^{\sharp} \right]_{ijkl} @@ -729,7 +729,7 @@ namespace Physics /** * Return the result of the pull back transformation on a rank-4 - * contravariant tensor, i.e. (in index notation) + * contravariant tensor, i.e. (in index notation): * @f[ * \textrm{det} \mathbf{F} \; \left[ * \chi^{-1}\left(\bullet\right)^{\sharp} \right]_{IJKL} @@ -750,7 +750,7 @@ namespace Physics /** * Return the result of the pull back transformation on a rank-4 - * contravariant symmetric tensor, i.e. (in index notation) + * contravariant symmetric tensor, i.e. (in index notation): * @f[ * \textrm{det} \mathbf{F} \; \left[ * \chi^{-1}\left(\bullet\right)^{\sharp} \right]_{IJKL} @@ -805,138 +805,98 @@ namespace Physics const Tensor<2, dim, Number> &F); //@} - } // namespace Transformations -} // namespace Physics - - -#ifndef DOXYGEN - -// ------------------------- inline functions ------------------------ + /** + * @name Basis transformations + */ + //@{ -namespace internal -{ - namespace Physics - { + /** + * Return a vector with a changed basis, i.e. + * @f[ + * \mathbf{V}^{\prime} \dealcoloneq \mathbf{B} \cdot \mathbf{V} + * @f] + * + * @param[in] V The vector to be transformed $\mathbf{V}$ + * @param[in] B The transformation matrix $\mathbf{B}$ + * @return $\mathbf{V}^{\prime}$ + */ template - inline Tensor<1, dim, Number> - transformation_contraction(const Tensor<1, dim, Number> &V, - const Tensor<2, dim, Number> &F) - { - return contract<1, 0>(F, V); - } - - + Tensor<1, dim, Number> + basis_transformation(const Tensor<1, dim, Number> &V, + const Tensor<2, dim, Number> &B); + /** + * Return a rank-2 tensor with a changed basis, i.e. + * @f[ + * \mathbf{T}^{\prime} \dealcoloneq \mathbf{B} \cdot \mathbf{T} \cdot + * \mathbf{B}^{T} + * @f] + * + * @param[in] T The tensor to be transformed $\mathbf{T}$ + * @param[in] B The transformation matrix $\mathbf{B}$ + * @return $\mathbf{T}^{\prime}$ + */ template - inline Tensor<2, dim, Number> - transformation_contraction(const Tensor<2, dim, Number> &T, - const Tensor<2, dim, Number> &F) - { - return contract<1, 0>(F, contract<1, 1>(T, F)); - } - - + Tensor<2, dim, Number> + basis_transformation(const Tensor<2, dim, Number> &T, + const Tensor<2, dim, Number> &B); + /** + * Return a symmetric rank-2 tensor with a changed basis, i.e. + * @f[ + * \mathbf{T}^{\prime} \dealcoloneq \mathbf{B} \cdot \mathbf{T} \cdot + * \mathbf{B}^{T} + * @f] + * + * @param[in] T The tensor to be transformed $\mathbf{T}$ + * @param[in] B The transformation matrix $\mathbf{B}$ + * @return $\mathbf{T}^{\prime}$ + */ template - inline dealii::SymmetricTensor<2, dim, Number> - transformation_contraction(const dealii::SymmetricTensor<2, dim, Number> &T, - const Tensor<2, dim, Number> & F) - { - Tensor<2, dim, Number> tmp_1; - for (unsigned int i = 0; i < dim; ++i) - for (unsigned int J = 0; J < dim; ++J) - // Loop over I but complex.h defines a macro I, so use I_ instead - for (unsigned int I_ = 0; I_ < dim; ++I_) - tmp_1[i][J] += F[i][I_] * T[I_][J]; - - dealii::SymmetricTensor<2, dim, Number> out; - for (unsigned int i = 0; i < dim; ++i) - for (unsigned int j = i; j < dim; ++j) - for (unsigned int J = 0; J < dim; ++J) - out[i][j] += F[j][J] * tmp_1[i][J]; - - return out; - } - - + SymmetricTensor<2, dim, Number> + basis_transformation(const SymmetricTensor<2, dim, Number> &T, + const Tensor<2, dim, Number> & B); + /** + * Return a rank-4 tensor with a changed basis, i.e. (in index notation): + * @f[ + * H_{ijkl}^{\prime} \dealcoloneq B_{iI} B_{jJ} H_{IJKL} B_{kK} B_{lL} + * @f] + * + * @param[in] H The tensor to be transformed $\mathbf{T}$ + * @param[in] B The transformation matrix $\mathbf{B}$ + * @return $\mathbf{H}^{\prime}$ + */ template - inline Tensor<4, dim, Number> - transformation_contraction(const Tensor<4, dim, Number> &H, - const Tensor<2, dim, Number> &F) - { - // This contraction order and indexing might look a bit dubious, so a - // quick explanation as to what's going on is probably in order: - // - // When the contract() function operates on the inner indices, the - // result has the inner index and outer index transposed, i.e. - // contract<2,1>(H,F) implies - // T_{IJLk} = (H_{IJMN} F_{mM}) \delta_{mL} \delta_{Nk} - // rather than T_{IJkL} (the desired result). - // So, in effect, contraction of the 3rd (inner) index with F as the - // second argument results in its transposition with respect to its - // adjacent neighbor. This is due to the position of the argument F, - // leading to the free index being on the right hand side of the result. - // However, given that we can do two transformations from the LHS of H - // and two from the right we can undo the otherwise erroneous - // swapping of the outer indices upon application of the second - // sets of contractions. - // - // Note: Its significantly quicker (in 3d) to push forward - // each index individually - return contract<1, 1>( - F, contract<1, 1>(F, contract<2, 1>(contract<2, 1>(H, F), F))); - } + Tensor<4, dim, Number> + basis_transformation(const Tensor<4, dim, Number> &H, + const Tensor<2, dim, Number> &B); + /** + * Return a symmetric rank-4 tensor with a changed basis, i.e. (in index + * notation): + * @f[ + * H_{ijkl}^{\prime} \dealcoloneq B_{iI} B_{jJ} H_{IJKL} B_{kK} B_{lL} + * @f] + * + * @param[in] H The tensor to be transformed $\mathbf{T}$ + * @param[in] B The transformation matrix $\mathbf{B}$ + * @return $\mathbf{H}^{\prime}$ + */ + template + SymmetricTensor<4, dim, Number> + basis_transformation(const SymmetricTensor<4, dim, Number> &H, + const Tensor<2, dim, Number> & B); + //@} - template - inline dealii::SymmetricTensor<4, dim, Number> - transformation_contraction(const dealii::SymmetricTensor<4, dim, Number> &H, - const Tensor<2, dim, Number> & F) - { - // The first and last transformation operations respectively - // break and recover the symmetry properties of the tensors. - // We also want to perform a minimal number of operations here - // and avoid some complications related to the transposition of - // tensor indices when contracting inner indices using the contract() - // function. (For an explanation of the contraction operations, - // please see the note in the equivalent function for standard - // Tensors.) So what we'll do here is manually perform the first - // and last contractions that break/recover the tensor symmetries - // on the inner indices, and use the contract() function only on - // the outer indices. - // - // Note: Its significantly quicker (in 3d) to push forward - // each index individually - - // Push forward (inner) index 1 - Tensor<4, dim, Number> tmp; - // Loop over I but complex.h defines a macro I, so use I_ instead - for (unsigned int I_ = 0; I_ < dim; ++I_) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int K = 0; K < dim; ++K) - for (unsigned int L = 0; L < dim; ++L) - for (unsigned int J = 0; J < dim; ++J) - tmp[I_][j][K][L] += F[j][J] * H[I_][J][K][L]; + } // namespace Transformations +} // namespace Physics - // Push forward (outer) indices 0 and 3 - tmp = contract<1, 0>(F, contract<3, 1>(tmp, F)); - // Push forward (inner) index 2 - dealii::SymmetricTensor<4, dim, Number> out; - for (unsigned int i = 0; i < dim; ++i) - for (unsigned int j = i; j < dim; ++j) - for (unsigned int k = 0; k < dim; ++k) - for (unsigned int l = k; l < dim; ++l) - for (unsigned int K = 0; K < dim; ++K) - out[i][j][k][l] += F[k][K] * tmp[i][j][K][l]; - return out; - } - } // namespace Physics -} // namespace internal +#ifndef DOXYGEN @@ -982,7 +942,7 @@ Physics::Transformations::Contravariant::push_forward( const Tensor<1, dim, Number> &V, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(V, F); + return Physics::Transformations::basis_transformation(V, F); } @@ -993,7 +953,7 @@ Physics::Transformations::Contravariant::push_forward( const Tensor<2, dim, Number> &T, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(T, F); + return Physics::Transformations::basis_transformation(T, F); } @@ -1004,7 +964,7 @@ Physics::Transformations::Contravariant::push_forward( const SymmetricTensor<2, dim, Number> &T, const Tensor<2, dim, Number> & F) { - return internal::Physics::transformation_contraction(T, F); + return Physics::Transformations::basis_transformation(T, F); } @@ -1015,7 +975,7 @@ Physics::Transformations::Contravariant::push_forward( const Tensor<4, dim, Number> &H, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(H, F); + return Physics::Transformations::basis_transformation(H, F); } @@ -1026,7 +986,7 @@ Physics::Transformations::Contravariant::push_forward( const SymmetricTensor<4, dim, Number> &H, const Tensor<2, dim, Number> & F) { - return internal::Physics::transformation_contraction(H, F); + return Physics::Transformations::basis_transformation(H, F); } @@ -1037,7 +997,7 @@ Physics::Transformations::Contravariant::pull_back( const Tensor<1, dim, Number> &v, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(v, invert(F)); + return Physics::Transformations::basis_transformation(v, invert(F)); } @@ -1048,7 +1008,7 @@ Physics::Transformations::Contravariant::pull_back( const Tensor<2, dim, Number> &t, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(t, invert(F)); + return Physics::Transformations::basis_transformation(t, invert(F)); } @@ -1059,7 +1019,7 @@ Physics::Transformations::Contravariant::pull_back( const SymmetricTensor<2, dim, Number> &t, const Tensor<2, dim, Number> & F) { - return internal::Physics::transformation_contraction(t, invert(F)); + return Physics::Transformations::basis_transformation(t, invert(F)); } @@ -1070,7 +1030,7 @@ Physics::Transformations::Contravariant::pull_back( const Tensor<4, dim, Number> &h, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(h, invert(F)); + return Physics::Transformations::basis_transformation(h, invert(F)); } @@ -1081,7 +1041,7 @@ Physics::Transformations::Contravariant::pull_back( const SymmetricTensor<4, dim, Number> &h, const Tensor<2, dim, Number> & F) { - return internal::Physics::transformation_contraction(h, invert(F)); + return Physics::Transformations::basis_transformation(h, invert(F)); } @@ -1092,7 +1052,8 @@ Physics::Transformations::Covariant::push_forward( const Tensor<1, dim, Number> &V, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(V, transpose(invert(F))); + return Physics::Transformations::basis_transformation(V, + transpose(invert(F))); } @@ -1103,7 +1064,8 @@ Physics::Transformations::Covariant::push_forward( const Tensor<2, dim, Number> &T, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(T, transpose(invert(F))); + return Physics::Transformations::basis_transformation(T, + transpose(invert(F))); } @@ -1114,7 +1076,8 @@ Physics::Transformations::Covariant::push_forward( const SymmetricTensor<2, dim, Number> &T, const Tensor<2, dim, Number> & F) { - return internal::Physics::transformation_contraction(T, transpose(invert(F))); + return Physics::Transformations::basis_transformation(T, + transpose(invert(F))); } @@ -1125,7 +1088,8 @@ Physics::Transformations::Covariant::push_forward( const Tensor<4, dim, Number> &H, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(H, transpose(invert(F))); + return Physics::Transformations::basis_transformation(H, + transpose(invert(F))); } @@ -1136,7 +1100,8 @@ Physics::Transformations::Covariant::push_forward( const SymmetricTensor<4, dim, Number> &H, const Tensor<2, dim, Number> & F) { - return internal::Physics::transformation_contraction(H, transpose(invert(F))); + return Physics::Transformations::basis_transformation(H, + transpose(invert(F))); } @@ -1146,7 +1111,7 @@ inline Tensor<1, dim, Number> Physics::Transformations::Covariant::pull_back(const Tensor<1, dim, Number> &v, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(v, transpose(F)); + return Physics::Transformations::basis_transformation(v, transpose(F)); } @@ -1156,7 +1121,7 @@ inline Tensor<2, dim, Number> Physics::Transformations::Covariant::pull_back(const Tensor<2, dim, Number> &t, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(t, transpose(F)); + return Physics::Transformations::basis_transformation(t, transpose(F)); } @@ -1167,7 +1132,7 @@ Physics::Transformations::Covariant::pull_back( const SymmetricTensor<2, dim, Number> &t, const Tensor<2, dim, Number> & F) { - return internal::Physics::transformation_contraction(t, transpose(F)); + return Physics::Transformations::basis_transformation(t, transpose(F)); } @@ -1177,7 +1142,7 @@ inline Tensor<4, dim, Number> Physics::Transformations::Covariant::pull_back(const Tensor<4, dim, Number> &h, const Tensor<2, dim, Number> &F) { - return internal::Physics::transformation_contraction(h, transpose(F)); + return Physics::Transformations::basis_transformation(h, transpose(F)); } @@ -1188,7 +1153,7 @@ Physics::Transformations::Covariant::pull_back( const SymmetricTensor<4, dim, Number> &h, const Tensor<2, dim, Number> & F) { - return internal::Physics::transformation_contraction(h, transpose(F)); + return Physics::Transformations::basis_transformation(h, transpose(F)); } @@ -1305,6 +1270,127 @@ Physics::Transformations::nansons_formula(const Tensor<1, dim, Number> &N, return cofactor(F) * N; } + +template +inline Tensor<1, dim, Number> +Physics::Transformations::basis_transformation(const Tensor<1, dim, Number> &V, + const Tensor<2, dim, Number> &B) +{ + return contract<1, 0>(B, V); +} + + + +template +inline Tensor<2, dim, Number> +Physics::Transformations::basis_transformation(const Tensor<2, dim, Number> &T, + const Tensor<2, dim, Number> &B) +{ + return contract<1, 0>(B, contract<1, 1>(T, B)); +} + + + +template +inline SymmetricTensor<2, dim, Number> +Physics::Transformations::basis_transformation( + const SymmetricTensor<2, dim, Number> &T, + const Tensor<2, dim, Number> & B) +{ + Tensor<2, dim, Number> tmp_1; + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int J = 0; J < dim; ++J) + // Loop over I but complex.h defines a macro I, so use I_ instead + for (unsigned int I_ = 0; I_ < dim; ++I_) + tmp_1[i][J] += B[i][I_] * T[I_][J]; + + SymmetricTensor<2, dim, Number> out; + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = i; j < dim; ++j) + for (unsigned int J = 0; J < dim; ++J) + out[i][j] += B[j][J] * tmp_1[i][J]; + + return out; +} + + + +template +inline Tensor<4, dim, Number> +Physics::Transformations::basis_transformation(const Tensor<4, dim, Number> &H, + const Tensor<2, dim, Number> &B) +{ + // This contraction order and indexing might look a bit dubious, so a + // quick explanation as to what's going on is probably in order: + // + // When the contract() function operates on the inner indices, the + // result has the inner index and outer index transposed, i.e. + // contract<2,1>(H,F) implies + // T_{IJLk} = (H_{IJMN} F_{mM}) \delta_{mL} \delta_{Nk} + // rather than T_{IJkL} (the desired result). + // So, in effect, contraction of the 3rd (inner) index with F as the + // second argument results in its transposition with respect to its + // adjacent neighbor. This is due to the position of the argument F, + // leading to the free index being on the right hand side of the result. + // However, given that we can do two transformations from the LHS of H + // and two from the right we can undo the otherwise erroneous + // swapping of the outer indices upon application of the second + // sets of contractions. + // + // Note: Its significantly quicker (in 3d) to push forward + // each index individually + return contract<1, 1>( + B, contract<1, 1>(B, contract<2, 1>(contract<2, 1>(H, B), B))); +} + + + +template +inline SymmetricTensor<4, dim, Number> +Physics::Transformations::basis_transformation( + const SymmetricTensor<4, dim, Number> &H, + const Tensor<2, dim, Number> & B) +{ + // The first and last transformation operations respectively + // break and recover the symmetry properties of the tensors. + // We also want to perform a minimal number of operations here + // and avoid some complications related to the transposition of + // tensor indices when contracting inner indices using the contract() + // function. (For an explanation of the contraction operations, + // please see the note in the equivalent function for standard + // Tensors.) So what we'll do here is manually perform the first + // and last contractions that break/recover the tensor symmetries + // on the inner indices, and use the contract() function only on + // the outer indices. + // + // Note: Its significantly quicker (in 3d) to push forward + // each index individually + + // Push forward (inner) index 1 + Tensor<4, dim, Number> tmp; + // Loop over I but complex.h defines a macro I, so use I_ instead + for (unsigned int I_ = 0; I_ < dim; ++I_) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int K = 0; K < dim; ++K) + for (unsigned int L = 0; L < dim; ++L) + for (unsigned int J = 0; J < dim; ++J) + tmp[I_][j][K][L] += B[j][J] * H[I_][J][K][L]; + + // Push forward (outer) indices 0 and 3 + tmp = contract<1, 0>(B, contract<3, 1>(tmp, B)); + + // Push forward (inner) index 2 + SymmetricTensor<4, dim, Number> out; + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = i; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = k; l < dim; ++l) + for (unsigned int K = 0; K < dim; ++K) + out[i][j][k][l] += B[k][K] * tmp[i][j][K][l]; + + return out; +} + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE