From c5afea76b13409eab09655fcc456c4635c2c1c67 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 11 Sep 2015 14:46:07 -0500 Subject: [PATCH] Implement TensorAccessors::contract3 --- include/deal.II/base/tensor_accessors.h | 107 ++++++++++++++++++++++++ 1 file changed, 107 insertions(+) diff --git a/include/deal.II/base/tensor_accessors.h b/include/deal.II/base/tensor_accessors.h index 02c6ab5589..2c0f4750d7 100644 --- a/include/deal.II/base/tensor_accessors.h +++ b/include/deal.II/base/tensor_accessors.h @@ -74,6 +74,7 @@ namespace TensorAccessors template class ReorderedIndexView; template struct ExtractHelper; template class Contract; + template class Contract3; } @@ -281,6 +282,41 @@ namespace TensorAccessors } + /** + * Full contraction of three tensorial objects: + * + * @f[ + * \sum_{i_1,..,i_{r1},j_1,..,j_{r2}} + * \text{left}_{i_1,..,i_{r1}} + * \text{middle}_{i_1,..,i_{r1},j_1,..,j_{r2}} + * \text{right}_{j_1,..,j_{r2}} + * @f] + * + * Calling this function is equivalent of writing the following low level + * code: + * @code + * T1 result = T1(); + * for(unsigned int i_0 = 0; i_0 < dim; ++i_0) + * ... + * for(unsigned int i_ = 0; i_ < dim; ++i_) + * for(unsigned int j_0 = 0; j_0 < dim; ++j_0) + * ... + * for(unsigned int j_ = 0; j_ < dim; ++j_) + * result += left[i_0]..[i_] * middle[i_0]..[i_][j_0]..[j_] * right[j_0]..[j_]; + * @endcode + * + * @note The Types @p T2, @p T3, and @p T4 must have + * rank rank_1, rank_1 + rank_2, and rank_3, respectively. @p T1 + * must be a scalar type. + */ + template + T1 contract3(const T2 &left, const T3 &middle, const T4 &right) + { + return internal::Contract3::template contract3 + (left, middle, right); + } + + namespace internal { // ------------------------------------------------------------------------- @@ -648,6 +684,77 @@ namespace TensorAccessors } }; + + // ------------------------------------------------------------------------- + // Implemenation of helper classes for contract3 + // ------------------------------------------------------------------------- + + // Fully contract three tensorial objects + // + // As long as rank_1 > 0, recurse over left and middle: + // + // for(unsigned int i_0; i_0 < dim; ++i_0) + // ... + // for(i_; i_ < dim; ++i_) + // [...] + // left[i_0]..[i_] ... middle[i_0]..[i_] ... right + + template + class Contract3 + { + public: + template + static inline + T1 contract3(const T2 &left, const T3 &middle, const T4 &right) + { + T1 result = T1(); + for (unsigned int i = 0; i < dim; ++i) + result += Contract3::template contract3(left[i], middle[i], right); + return result; + } + }; + + // If rank_1 ==0, continue to recurse over middle and right: + // + // for(unsigned int i_0; i_0 < dim; ++i_0) + // ... + // for(i_; i_ < dim; ++i_) + // for(unsigned int j_0; j_0 < dim; ++j_0) + // ... + // for(j_; j_ < dim; ++j_) + // [...] + // left[i_0]..[i_] ... middle[i_0]..[i_][j_0]..[j_] ... right[j_0]..[j_] + + template + class Contract3<0, rank_2, dim> + { + public: + template + static inline + T1 contract3(const T2 &left, const T3 &middle, const T4 &right) + { + T1 result = T1(); + for (unsigned int i = 0; i < dim; ++i) + result += Contract3<0, rank_2 - 1, dim>::template contract3(left, middle[i], right[i]); + return result; + } + }; + + // Contraction of three tensorial objects of rank 0 is just a scalar + // multiplication. + + template + class Contract3<0, 0, dim> + { + public: + template + static inline + T1 contract3(const T2 &left, const T3 &middle, const T4 &right) + { + return left * middle * right; + } + }; + // ------------------------------------------------------------------------- } /* namespace internal */ -- 2.39.5