]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add contract3() function with SymmetricTensor arguments.
authorJean-Paul Pelteret <jppelteret@gmail.com>
Tue, 10 Jan 2017 15:10:26 +0000 (16:10 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 12 Jan 2017 15:17:58 +0000 (16:17 +0100)
fixes #2691

doc/news/changes/minor/20170110Jean-PaulPelteretMatthiasMaier [new file with mode: 0644]
include/deal.II/base/symmetric_tensor.h
include/deal.II/base/tensor.h
tests/base/symmetric_tensor_37.cc [new file with mode: 0644]
tests/base/symmetric_tensor_37.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170110Jean-PaulPelteretMatthiasMaier b/doc/news/changes/minor/20170110Jean-PaulPelteretMatthiasMaier
new file mode 100644 (file)
index 0000000..1d7bf7f
--- /dev/null
@@ -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)
index ae55560d157479b3bc436f4ef7a4dd79b654908d..86fd420cfdd3a948c1a00d3745da2f6024403839 100644 (file)
@@ -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));
+    }
   }
 }
 
index 340e5bb2011e99de3b542aeb2f4b86da9ba794db..da9fadd3d89fb45a42fbb6e9c821eb5732cb0aef 100644 (file)
@@ -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 (file)
index 0000000..8776671
--- /dev/null
@@ -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 (file)
index 0000000..f636124
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::All OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.