From c42a780991bcb7dbaaea672b6ddccf0c5cefa0d4 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret 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. +
+(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 operator [] (const unsigned int i); + /** + * Index operator. + */ + Accessor 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 + Accessor + Accessor::operator[] (const unsigned int i) const + { + return Accessor (tensor, + merge (previous_indices, i, rank-P)); + } + + + template Accessor:: Accessor (tensor_type &tensor, @@ -894,6 +914,14 @@ namespace internal { return tensor(merge (previous_indices, i, rank-1)); } + + + template + typename Accessor::reference + Accessor::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 &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 class TensorT1, + template class TensorT2, + template class TensorT3, + int rank_1, int rank_2, int dim, typename T1, typename T2, typename T3> typename ProductType::type>::type -contract3 (const Tensor &left, - const Tensor &middle, - const Tensor &right) +contract3 (const TensorT1 &left, + const TensorT2 &middle, + const TensorT3 &right) { typedef typename ProductType::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 +#include +#include + +#include +#include + +// Although not strictly required, in this test it is expected that T1 or T3 +// (or both) are a symmetric tensor +template class T1, + template class T3> +void test_symm_tensor_contract_3(const T1 &l, + const Tensor &m, + const T3 &r) +{ + const double res1 = contract3(l,m,r); + const double res2 = contract3(static_cast >(l), + m, + static_cast >(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