From db9574b049ff6a7e3dcd1ef913a61364b37635af Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Tue, 17 Oct 2017 13:28:41 +0200 Subject: [PATCH] Fix Physics::internal::transformation_contraction for standard Tensors. A robust unit test has been attached free of charge. --- .../changes/minor/20171017Jean-PaulPelteret | 5 + include/deal.II/physics/transformations.h | 66 +++-- .../transformations-push_forward_01.cc | 248 ++++++++++++++++++ .../transformations-push_forward_01.output | 2 + 4 files changed, 297 insertions(+), 24 deletions(-) create mode 100644 doc/news/changes/minor/20171017Jean-PaulPelteret create mode 100644 tests/physics/transformations-push_forward_01.cc create mode 100644 tests/physics/transformations-push_forward_01.output diff --git a/doc/news/changes/minor/20171017Jean-PaulPelteret b/doc/news/changes/minor/20171017Jean-PaulPelteret new file mode 100644 index 0000000000..22554691f6 --- /dev/null +++ b/doc/news/changes/minor/20171017Jean-PaulPelteret @@ -0,0 +1,5 @@ +Fixed: The transformation operation for rank-2 and rank-4 (non-symmetric) tensors +in the Physics::Transformations was incorrectly implemented. +This has now been fixed. +
+(Jean-Paul Pelteret, 2017/10/17) diff --git a/include/deal.II/physics/transformations.h b/include/deal.II/physics/transformations.h index 88de5ef148..b5c7a0a065 100644 --- a/include/deal.II/physics/transformations.h +++ b/include/deal.II/physics/transformations.h @@ -818,7 +818,7 @@ namespace internal transformation_contraction (const Tensor<2,dim,Number> &T, const Tensor<2,dim,Number> &F) { - return contract<1,1>(F,contract<1,0>(F, T)); + return contract<1,0>(F,contract<1,1>(T, F)); } @@ -852,9 +852,26 @@ namespace internal transformation_contraction (const Tensor<4,dim,Number> &H, const Tensor<2,dim,Number> &F) { - // Its significantly quicker (in 3d) to push forward + // 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,3>(F,contract<1,2>(F,contract<1,1>(F,contract<1,0>(F, H)))); + return contract<1,1>(F,contract<1,1>(F,contract<2,1>(contract<2,1>(H,F), F))); } @@ -865,40 +882,41 @@ namespace internal transformation_contraction (const dealii::SymmetricTensor<4,dim,Number> &H, const Tensor<2,dim,Number> &F) { - // Its significantly quicker (in 3d) to push forward + // 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 - Tensor<4,dim,Number> tmp_1; - for (unsigned int i=0; i tmp_2; - for (unsigned int i=0; i tmp; + for (unsigned int I=0; I(F, contract<3,1>(tmp,F)); + // Push forward (inner) index 2 dealii::SymmetricTensor<4,dim,Number> out; for (unsigned int i=0; i +#include +#include + +using namespace dealii; + +template +void test_tensor (const Tensor<2,dim> &F) +{ + // Rank-1 Tensors + { + Tensor<1,dim> T; + unsigned int c=1; + for (unsigned int i=0; i T_calc; + for (unsigned int i=0; i T_trans = Physics::Transformations::Contravariant::push_forward(T,F); + Assert((T_calc - T_trans).norm() < 1e-9, ExcMessage("Rank 1 tensor: Contraction using push_forward() function is incorrect.")); + } + + // Rank-2 Tensors + { + Tensor<2,dim> T; + unsigned int c=1; + for (unsigned int i=0; i tmp0; + for (unsigned int i=0; i T_calc; + for (unsigned int i=0; i T_trans = Physics::Transformations::Contravariant::push_forward(T,F); + Assert((T_calc - T_trans).norm() < 1e-9, ExcMessage("Rank 2 tensor: Contraction using push_forward() function is incorrect.")); + } + + // Rank-4 Tensors + { + Tensor<4,dim> T; + unsigned int c=1; + for (unsigned int i=0; i tmp0; + for (unsigned int i=0; i tmp1; + for (unsigned int i=0; i tmp2; + for (unsigned int i=0; i T_calc; + for (unsigned int i=0; i T_trans = Physics::Transformations::Contravariant::push_forward(T,F); + Assert((T_calc - T_trans).norm() < 1e-9, ExcMessage("Rank 4 tensor: Contraction using push_forward() function is incorrect.")); + } +} + +template +void test_symmetric_tensor (const Tensor<2,dim> &F) +{ + // Rank-2 Symmetric tensors + { + SymmetricTensor<2,dim> T; + unsigned int c=1; + for (unsigned int i=0; i tmp0; + for (unsigned int i=0; i T_calc; + for (unsigned int i=0; i T_trans = Physics::Transformations::Contravariant::push_forward(T,F); + Assert((T_calc - T_trans).norm() < 1e-9, ExcMessage("Rank 2 symmetric tensor: Contraction using push_forward() function is incorrect.")); + } + + // Rank-4 Symmetric tensors + { + SymmetricTensor<4,dim> T; + unsigned int c=1; + for (unsigned int i=0; i tmp0; + for (unsigned int i=0; i tmp1; + for (unsigned int i=0; i tmp2; + for (unsigned int i=0; i T_calc; + for (unsigned int i=0; i T_trans = Physics::Transformations::Contravariant::push_forward(T,F); + + std::cout << "T_calc: " << T_calc << std::endl; + std::cout << "T_trans: " << T_trans << std::endl; + + Assert((T_calc - T_trans).norm() < 1e-9, ExcMessage("Rank 4 symmetric tensor: Contraction using push_forward() function is incorrect.")); + } +} + +template +void test () +{ + // Test with unit tensor + test_tensor(unit_symmetric_tensor()); + test_symmetric_tensor(unit_symmetric_tensor()); + + // Test with non-trivial tensor + Tensor<2,dim> F = unit_symmetric_tensor(); + double c=0.1; + for (unsigned int i=0; i(F); + test_symmetric_tensor(F); +} + +int main (int argc, char *argv[]) +{ + initlog(); + + test<2>(); + test<3>(); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/physics/transformations-push_forward_01.output b/tests/physics/transformations-push_forward_01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/physics/transformations-push_forward_01.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5