--- /dev/null
+Improved: The contract3() function has been extended accommodate both Tensor
+and SymmetricTensors arguments.
+<br>
+(Jean-Paul Pelteret, Matthias Maier, 2017/01/10)
*/
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.
*/
reference operator [] (const unsigned int);
+ /**
+ * Index operator.
+ */
+ reference operator [] (const unsigned int) const;
+
private:
/**
* Store the data given to the constructor.
+ 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,
{
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));
+ }
}
}
* \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;
--- /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 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;
+}
--- /dev/null
+
+DEAL::All OK