From 39a32b594609164449d6835e127400db8167de8c Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Wed, 16 Sep 2015 20:52:19 -0500 Subject: [PATCH] Implement a generic contract function for contraction of one or two indices --- include/deal.II/base/tensor.h | 143 +++++++++++++++++++++++++++++- tests/base/symmetric_tensor_06.cc | 2 +- 2 files changed, 143 insertions(+), 2 deletions(-) diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index 031f1cd5e3..bb5e1edf28 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -1487,6 +1487,147 @@ operator * (const Tensor &src1, } +/** + * Generic contraction of a pair of indices of two tensors of arbitrary + * rank: Return a tensor of rank $(\text{rank_1} + \text{rank_2} - 2)$ that + * is the contraction of index @p index_1 of a tensor @p src1 of rank + * @p rank_1 with the index @p index_2 of a tensor @p src2 of rank @p rank_2: + * @f[ + * \text{result}_{i_1,..,i_{r1},j_1,..,j_{r2}} + * = \sum_{k} + * \text{left}_{i_1,..,k,..,i_{r1}} + * \text{right}_{j_1,..,k,..,j_{r2}} + * @f] + * + * If for example the first index (index_1==0) of a tensor + * t1 shall be contracted with the third index + * (index_2==2) of a tensor t2, the invocation of + * this function is + * @code + * contract<0, 2>(t1, t2); + * @endcode + * + * @note The position of the index is counted from 0, i.e., + * $0\le\text{index_i}<\text{range_i}$. + * + * @note In case the contraction yields tensor of rank 0 the scalar + * number is returned as an unwrapped number type. + * + * @relates Tensor + * @author Matthias Maier, 2015 + */ +template +inline +typename Tensor::type>::tensor_type +contract (const Tensor &src1, + const Tensor &src2) +{ +#ifdef DEAL_II_WITH_CXX11 + static_assert(0 <= index_1 && index_1 < rank_1, + "The specified index_1 must lie within the range [0,rank_1)"); + static_assert(0 <= index_2 && index_2 < rank_2, + "The specified index_2 must lie within the range [0,rank_2)"); +#endif + using namespace TensorAccessors; + using namespace TensorAccessors::internal; + + // Reorder index_1 to the end of src1: + ReorderedIndexView > + reord_01 = reordered_index_view(src1); + + // Reorder index_2 to the end of src2: + ReorderedIndexView > + reord_02 = reordered_index_view(src2); + + typename Tensor::type>::tensor_type + result; + TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02); + return result; +} + + +/** + * Generic contraction of two pairs of indices of two tensors of + * arbitrary rank: Return a tensor of rank + * $(\text{rank_1} + \text{rank_2} - 4)$ that is the contraction of index + * @p index_1 with index @p index_2, and index @p index_3 with index + * @p index_4 of a tensor @p src1 of rank @p rank_1 and a tensor @p src2 of + * rank @p rank_2: + * @f[ + * \text{result}_{i_1,..,i_{r1},j_1,..,j_{r2}} + * = \sum_{k, l} + * \text{left}_{i_1,..,k,..,l,..,i_{r1}} + * \text{right}_{j_1,..,k,..,l..,j_{r2}} + * @f] + * + * If for example the first index (index_1==0) shall be + * contracted with the third index (index_2==2), and the + * second index (index_3==1) with the first index + * (index_4==0) the invocation of this function is + * this function is + * @code + * contract<0, 2, 1, 0>(t1, t2); + * @endcode + * + * @note The position of the index is counted from 0, i.e., + * $0\le\text{index_i}<\text{range_i}$. + * + * @note In case the contraction yields tensor of rank 0 the scalar + * number is returned as an unwrapped number type. + * + * @relates Tensor + * @author Matthias Maier, 2015 + */ +template +inline +typename Tensor::type>::tensor_type +contract (const Tensor &src1, + const Tensor &src2) +{ +#ifdef DEAL_II_WITH_CXX11 + static_assert(0 <= index_1 && index_1 < rank_1, + "The specified index_1 must lie within the range [0,rank_1)"); + static_assert(0 <= index_3 && index_3 < rank_1, + "The specified index_3 must lie within the range [0,rank_1)"); + static_assert(0 <= index_2 && index_2 < rank_2, + "The specified index_2 must lie within the range [0,rank_2)"); + static_assert(0 <= index_4 && index_4 < rank_2, + "The specified index_4 must lie within the range [0,rank_2)"); +#endif + using namespace TensorAccessors; + using namespace TensorAccessors::internal; + + // Reorder index_1 to the end of src1: + ReorderedIndexView > + reord_1 = TensorAccessors::reordered_index_view(src1); + + // Reorder index_2 to the end of src2: + ReorderedIndexView > + reord_2 = TensorAccessors::reordered_index_view(src2); + + // Now, reorder index_3 to the end of src1. We have to make sure to + // preserve the orginial ordering: index_1 has been removed. If + // index_3 > index_1, we have to use (index_3 - 1) instead: + ReorderedIndexView<(index_3 < index_1 ? index_3 : index_3 - 1), rank_1, ReorderedIndexView > > + reord_3 = TensorAccessors::reordered_index_view(reord_1); + + // Now, reorder index_4 to the end of src2. We have to make sure to + // preserve the orginial ordering: index_2 has been removed. If + // index_4 > index_2, we have to use (index_4 - 1) instead: + ReorderedIndexView<(index_4 < index_2 ? index_4 : index_4 - 1), rank_2, ReorderedIndexView > > + reord_4 = TensorAccessors::reordered_index_view(reord_2); + + typename Tensor::type>::tensor_type + result; + TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4); + return result; +} + + /** * The scalar product, or (generalized) Frobenius inner product of two * tensors of equal rank: Return a scalar number that is the result of a @@ -1844,7 +1985,7 @@ linfty_norm (const Tensor<2,dim,Number> &t) DEAL_II_NAMESPACE_CLOSE // include deprecated non-member functions operating on Tensor -// #include +#include #endif diff --git a/tests/base/symmetric_tensor_06.cc b/tests/base/symmetric_tensor_06.cc index 3494f2280a..2401e8c352 100644 --- a/tests/base/symmetric_tensor_06.cc +++ b/tests/base/symmetric_tensor_06.cc @@ -67,7 +67,7 @@ void test () } bs[i][j] = tmp_ij; } - double_contract (ba, ta, aa); + ba = contract<2,0,3,1>(ta, aa); for (unsigned int i=0; i