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));
}
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)));
}
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;
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
projects / dealii.git / commitdiff