From c42a780991bcb7dbaaea672b6ddccf0c5cefa0d4 Mon Sep 17 00:00:00 2001
From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Tue, 10 Jan 2017 16:10:26 +0100
Subject: [PATCH] Add contract3() function with SymmetricTensor arguments.

fixes #2691
---
 .../20170110Jean-PaulPelteretMatthiasMaier    |  4 +
 include/deal.II/base/symmetric_tensor.h       | 28 ++++++
 include/deal.II/base/tensor.h                 | 16 ++--
 tests/base/symmetric_tensor_37.cc             | 87 +++++++++++++++++++
 tests/base/symmetric_tensor_37.output         |  2 +
 5 files changed, 132 insertions(+), 5 deletions(-)
 create mode 100644 doc/news/changes/minor/20170110Jean-PaulPelteretMatthiasMaier
 create mode 100644 tests/base/symmetric_tensor_37.cc
 create mode 100644 tests/base/symmetric_tensor_37.output

diff --git a/doc/news/changes/minor/20170110Jean-PaulPelteretMatthiasMaier b/doc/news/changes/minor/20170110Jean-PaulPelteretMatthiasMaier
new file mode 100644
index 0000000000..1d7bf7fbf5
--- /dev/null
+++ b/doc/news/changes/minor/20170110Jean-PaulPelteretMatthiasMaier
@@ -0,0 +1,4 @@
+Improved: The contract3() function has been extended accommodate both Tensor
+and SymmetricTensors arguments.
+<br>
+(Jean-Paul Pelteret, Matthias Maier, 2017/01/10)
diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h
index ae55560d15..86fd420cfd 100644
--- a/include/deal.II/base/symmetric_tensor.h
+++ b/include/deal.II/base/symmetric_tensor.h
@@ -323,6 +323,11 @@ namespace internal
        */
       Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);
 
+      /**
+       * Index operator.
+       */
+      Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i) const;
+
     private:
       /**
        * Store the data given to the constructor.
@@ -406,6 +411,11 @@ namespace internal
        */
       reference operator [] (const unsigned int);
 
+      /**
+       * Index operator.
+       */
+      reference operator [] (const unsigned int) const;
+
     private:
       /**
        * Store the data given to the constructor.
@@ -867,6 +877,16 @@ namespace internal
 
 
 
+    template <int rank, int dim, bool constness, int P, typename Number>
+    Accessor<rank,dim,constness,P-1,Number>
+    Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i) const
+    {
+      return Accessor<rank,dim,constness,P-1,Number> (tensor,
+                                                      merge (previous_indices, i, rank-P));
+    }
+
+
+
     template <int rank, int dim, bool constness, typename Number>
     Accessor<rank,dim,constness,1,Number>::
     Accessor (tensor_type              &tensor,
@@ -894,6 +914,14 @@ namespace internal
     {
       return tensor(merge (previous_indices, i, rank-1));
     }
+
+
+    template <int rank, int dim, bool constness, typename Number>
+    typename Accessor<rank,dim,constness,1,Number>::reference
+    Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i) const
+    {
+      return tensor(merge (previous_indices, i, rank-1));
+    }
   }
 }
 
diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index 340e5bb201..da9fadd3d8 100644
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -1658,15 +1658,21 @@ scalar_product (const Tensor<rank, dim, Number> &left,
  *   \text{right}_{j_1,..,j_{r2}}
  * @f]
  *
+ * @note Each of the three input tensors can be either a Tensor or
+ * SymmetricTensor.
+ *
  * @relates Tensor
- * @author Matthias Maier, 2015
+ * @author Matthias Maier, 2015, Jean-Paul Pelteret 2017
  */
-template <int rank_1, int rank_2, int dim,
+template <template<int, int, typename> class TensorT1,
+          template<int, int, typename> class TensorT2,
+          template<int, int, typename> class TensorT3,
+          int rank_1, int rank_2, int dim,
           typename T1, typename T2, typename T3>
 typename ProductType<T1, typename ProductType<T2, T3>::type>::type
-contract3 (const Tensor<rank_1, dim, T1> &left,
-           const Tensor<rank_1 + rank_2, dim, T2> &middle,
-           const Tensor<rank_2, dim, T3> &right)
+contract3 (const TensorT1<rank_1, dim, T1>          &left,
+           const TensorT2<rank_1 + rank_2, dim, T2> &middle,
+           const TensorT3<rank_2, dim, T3>          &right)
 {
   typedef typename ProductType<T1, typename ProductType<T2, T3>::type>::type
   return_type;
diff --git a/tests/base/symmetric_tensor_37.cc b/tests/base/symmetric_tensor_37.cc
new file mode 100644
index 0000000000..87766712da
--- /dev/null
+++ b/tests/base/symmetric_tensor_37.cc
@@ -0,0 +1,87 @@
+// ---------------------------------------------------------------------
+//
+// 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 the contract3 function involving SymmetricTensors
+
+#include "../tests.h"
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/logstream.h>
+
+#include <fstream>
+#include <iomanip>
+
+// Although not strictly required, in this test it is expected that T1 or T3
+// (or both) are a symmetric tensor
+template<int rank1, int rank2, int rank3, int dim, typename number,
+         template<int,int,typename> class T1,
+         template<int,int,typename> class T3>
+void test_symm_tensor_contract_3(const T1<rank1,dim,number>     &l,
+                                 const Tensor<rank2,dim,number> &m,
+                                 const T3<rank3,dim,number>     &r)
+{
+  const double res1 = contract3(l,m,r);
+  const double res2 = contract3(static_cast<Tensor<rank1,dim> >(l),
+                                m,
+                                static_cast<Tensor<rank3,dim> >(r));
+  Assert(std::abs(res1 - res2) < 1e-12,
+         ExcMessage("Result from symmetric tensor contract3 is incorrect."));
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(5);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  const int dim = 3;
+  Tensor<1,dim> v1;
+  Tensor<2,dim> t1;
+  SymmetricTensor<2,dim> s1,s2;
+  Tensor<3,dim> T1;
+  Tensor<4,dim> H1;
+
+  for (unsigned int i=0; i<dim; ++i)
+    {
+      v1[i] = 1+i;
+      for (unsigned int j=0; j<dim; ++j)
+        {
+          t1[i][j] = 1+j+i*dim;
+          s1[i][j] = 1+j+i*dim;
+          s2[i][j] = 2+(j+i)*dim;
+          for (unsigned int k=0; k<dim; ++k)
+            {
+              T1[i][j][k] = 1+k+j*dim+i*dim*dim;
+              for (unsigned int l=0; l<dim; ++l)
+                H1[i][j][k][l] = 1+l+k*dim+j*dim*dim+i*dim*dim*dim;
+            }
+        }
+    }
+
+  // Vector + SymmetricTensor
+  test_symm_tensor_contract_3(v1,T1,s1);
+  test_symm_tensor_contract_3(s1,T1,v1);
+
+  // Tensor + SymmetricTensor
+  test_symm_tensor_contract_3(t1,H1,s1);
+  test_symm_tensor_contract_3(s1,H1,t1);
+
+  // SymmetricTensor + SymmetricTensor
+  test_symm_tensor_contract_3(s1,H1,s2);
+
+  deallog << "All OK" << std::endl;
+}
diff --git a/tests/base/symmetric_tensor_37.output b/tests/base/symmetric_tensor_37.output
new file mode 100644
index 0000000000..f636124a17
--- /dev/null
+++ b/tests/base/symmetric_tensor_37.output
@@ -0,0 +1,2 @@
+
+DEAL::All OK
-- 
2.39.5