From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Tue, 17 Oct 2017 11:28:41 +0000 (+0200)
Subject: Fix Physics::internal::transformation_contraction for standard Tensors.
X-Git-Tag: v9.0.0-rc1~939^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=db9574b049ff6a7e3dcd1ef913a61364b37635af;p=dealii.git

Fix Physics::internal::transformation_contraction for standard Tensors.

A robust unit test has been attached free of charge.
---

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