From db9574b049ff6a7e3dcd1ef913a61364b37635af Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret <jppelteret@gmail.com> 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. +<br> +(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<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 I=0; I<dim; ++I) - tmp_1[i][J][K][L] += F[i][I] * H[I][J][K][L]; - - Tensor<4,dim,Number> tmp_2; - for (unsigned int i=0; i<dim; ++i) + // Push forward (inner) index 1 + Tensor<4,dim,Number> tmp; + 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_2[i][j][K][L] += F[j][J] * tmp_1[i][J][K][L]; + tmp[I][j][K][L] += F[j][J] * H[I][J][K][L]; - tmp_1 = 0.0; - 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 K=0; K<dim; ++K) - tmp_1[i][j][k][L] += F[k][K] * tmp_2[i][j][K][L]; + // 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 L=0; L<dim; ++L) - out[i][j][k][l] += F[l][L] * tmp_1[i][j][k][L]; + for (unsigned int K=0; K<dim; ++K) + out[i][j][k][l] += F[k][K] * tmp[i][j][K][l]; return out; } diff --git a/tests/physics/transformations-push_forward_01.cc b/tests/physics/transformations-push_forward_01.cc new file mode 100644 index 0000000000..4113ef3055 --- /dev/null +++ b/tests/physics/transformations-push_forward_01.cc @@ -0,0 +1,248 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Check that the push-forward operation in the Physics::Transformations +// namespace works as expected for all tensor types. +// This gives particular scrutiny to the use of the contract() function +// of the tensor class. + +#include "../tests.h" +#include <deal.II/base/tensor.h> +#include <deal.II/base/symmetric_tensor.h> +#include <deal.II/physics/transformations.h> + +using namespace dealii; + +template <int dim> +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<dim; ++i) + T[i] = c++; + + // Hand calculation: Push forward index 0 + Tensor<1,dim> T_calc; + for (unsigned int i=0; i<dim; ++i) + for (unsigned int I=0; I<dim; ++I) + T_calc[i] += F[i][I]*T[I]; + + const Tensor<1,dim> 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<dim; ++i) + for (unsigned int j=0; j<dim; ++j) + T[i][j] = c++; + + // Hand calculation: Push forward index 0 + Tensor<2,dim> tmp0; + for (unsigned int i=0; i<dim; ++i) + for (unsigned int J=0; J<dim; ++J) + for (unsigned int I=0; I<dim; ++I) + tmp0[i][J] += F[i][I]*T[I][J]; + + // Hand calculation: Push forward index 1 + Tensor<2,dim> T_calc; + for (unsigned int i=0; i<dim; ++i) + for (unsigned int j=0; j<dim; ++j) + for (unsigned int J=0; J<dim; ++J) + T_calc[i][j] += F[j][J]*tmp0[i][J]; + + const Tensor<2,dim> 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<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) + T[i][j][k][l] = c++; + + // Hand calculation: Push forward index 0 + Tensor<4,dim> tmp0; + 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 I=0; I<dim; ++I) + tmp0[i][J][K][L] += F[i][I]*T[I][J][K][L]; + + // Hand calculation: Push forward index 1 + Tensor<4,dim> tmp1; + 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) + tmp1[i][j][K][L] += F[j][J]*tmp0[i][J][K][L]; + + // Hand calculation: Push forward index 2 + Tensor<4,dim> tmp2; + 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 K=0; K<dim; ++K) + tmp2[i][j][k][L] += F[k][K]*tmp1[i][j][K][L]; + + // Hand calculation: Push forward index 3 + Tensor<4,dim> T_calc; + 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 L=0; L<dim; ++L) + T_calc[i][j][k][l] += F[l][L]*tmp2[i][j][k][L]; + + const Tensor<4,dim> 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 <int dim> +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<dim; ++i) + for (unsigned int j=i; j<dim; ++j) // symmetry + T[i][j] = c++; + + // Hand calculation: Push forward index 0 + // Note: The symmetry of the initial tensor is lost after the + // transformation of the first index + Tensor<2,dim> tmp0; + for (unsigned int i=0; i<dim; ++i) + for (unsigned int J=0; J<dim; ++J) + for (unsigned int I=0; I<dim; ++I) + tmp0[i][J] += F[i][I]*T[I][J]; + + // Hand calculation: Push forward index 1 + // Note: At this point we recover the symmetry of the + // transformed tensor. + SymmetricTensor<2,dim> T_calc; + for (unsigned int i=0; i<dim; ++i) + for (unsigned int j=i; j<dim; ++j) // symmetry + for (unsigned int J=0; J<dim; ++J) + T_calc[i][j] += F[j][J]*tmp0[i][J]; + + const Tensor<2,dim> 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<dim; ++i) + for (unsigned int j=i; j<dim; ++j) // symmetry + for (unsigned int k=0; k<dim; ++k) + for (unsigned int l=k; l<dim; ++l) + T[i][j][k][l] = c++; + + // Hand calculation: Push forward index 0 + // Note: The symmetry of the initial tensor is lost after the + // transformation of the first index + Tensor<4,dim> tmp0; + 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 I=0; I<dim; ++I) + tmp0[i][J][K][L] += F[i][I]*T[I][J][K][L]; + + // Hand calculation: Push forward index 1 + Tensor<4,dim> tmp1; + 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) + tmp1[i][j][K][L] += F[j][J]*tmp0[i][J][K][L]; + + // Hand calculation: Push forward index 2 + Tensor<4,dim> tmp2; + 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 K=0; K<dim; ++K) + tmp2[i][j][k][L] += F[k][K]*tmp1[i][j][K][L]; + + // Hand calculation: Push forward index 3 + // Note: At this point we recover the symmetry of the + // transformed tensor. + SymmetricTensor<4,dim> T_calc; + for (unsigned int i=0; i<dim; ++i) + for (unsigned int j=i; j<dim; ++j) // symmetry + for (unsigned int k=0; k<dim; ++k) + for (unsigned int l=k; l<dim; ++l) // symmetry + for (unsigned int L=0; L<dim; ++L) + T_calc[i][j][k][l] += F[l][L]*tmp2[i][j][k][L]; + + const SymmetricTensor<4,dim> 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 <int dim> +void test () +{ + // Test with unit tensor + test_tensor<dim>(unit_symmetric_tensor<dim>()); + test_symmetric_tensor<dim>(unit_symmetric_tensor<dim>()); + + // Test with non-trivial tensor + Tensor<2,dim> F = unit_symmetric_tensor<dim>(); + double c=0.1; + for (unsigned int i=0; i<dim; ++i) + for (unsigned int j=0; j<dim; ++j) + { + F[i][j] += c; + c += 0.05; + } + test_tensor<dim>(F); + test_symmetric_tensor<dim>(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