From 6dbfa5db5b81520bbb501236624269875ff4aef9 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 11 Sep 2015 14:41:54 -0500 Subject: [PATCH] Implement TensorAccessors::reordered_index_view --- include/deal.II/base/tensor_accessors.h | 345 ++++++++++++++++++++++++ 1 file changed, 345 insertions(+) diff --git a/include/deal.II/base/tensor_accessors.h b/include/deal.II/base/tensor_accessors.h index 22e9fcf6f3..9d4e405a7f 100644 --- a/include/deal.II/base/tensor_accessors.h +++ b/include/deal.II/base/tensor_accessors.h @@ -68,8 +68,353 @@ DEAL_II_NAMESPACE_OPEN */ namespace TensorAccessors { + // forward declarations + namespace internal + { + template class ReorderedIndexView; + } + /** + * This class provides a local typedef @p value_type denoting the + * resulting type of an access with operator[](unsigned int). More + * precisely, @p value_type will be + * - T::value_type if T is a tensorial class providing a + * typedef value_type and does not have a const qualifier. + * - const T::value_type if T is a tensorial class + * providing a typedef value_type and does have a const + * qualifier. + * - const T::value_type if T is a tensorial class + * providing a typedef value_type and does have a const + * qualifier. + * - A if T is of array type A[...] + * - const A if T is of array type A[...] and + * does have a const qualifier. + */ + template + struct ValueType + { + typedef typename T::value_type value_type; + }; + + template + struct ValueType + { + typedef const typename T::value_type value_type; + }; + + template + struct ValueType + { + typedef T value_type; + }; + + template + struct ValueType + { + typedef const T value_type; + }; + + + /** + * This class provides a local typedef @p value_type that is equal to + * the typedef value_type after @p deref_steps + * recursive dereferenciations via ```operator[](unsigned int)```. + * Further, constness is preserved via the ValueType + * type trait, i.e., if T is const, ReturnType::value_type + * will also be const. + */ + template + struct ReturnType + { + typedef typename ReturnType::value_type>::value_type value_type; + }; + + template + struct ReturnType<0, T> + { + typedef T value_type; + }; + + + /** + * Provide a "tensorial view" to a reference @p t of a tensor object of + * rank @p rank in which the index @p index is shifted to the + * end. As an example consider a tensor of 5th order in dim=5 space + * dimensions that can be accessed through 5 recursive + * operator[]() invocations: + * @code + * Tensor<5, dim> tensor; + * tensor[0][1][2][3][4] = 42.; + * @endcode + * Index 1 (the 2nd index, count starts at 0) can now be shifted to the + * end via + * @code + * auto tensor_view = reordered_index_view<1, 5>(tensor); + * tensor_view[0][2][3][4][1] == 42.; // is true + * @endcode + * The usage of the dealii::Tensor type was solely for the sake of an + * example. The mechanism implemented by this function is available for + * fairly general tensorial types @p T. + * + * The purpose of this reordering facility is to be able to contract over + * an arbitrary index of two (ore more) tensors: + * - reorder the indices in mind to the end of the tensors + * - use belows contract function that contracts the _last_ elements of + * tensors. + * + * @note This function returns an internal class object consisting of an + * array subscript operator operator[](unsigned int) and a + * typedef value_type describing its return value. + * + * @tparam index The index to be shifted to the end. Indices are counted + * from 0, thus the valid range is $0\le\text{index}<\text{rank}$. + * @tparam rank Rank of the tensorial object @param t + * @tparam T A tensorial object of rank @p rank. @p T must + * provide a local typedef value_type and an index operator + * operator[]() that returns a (const or non-const) + * reference of value_type: + * @code + * class T + * { + * typedef ... value_type + * value_type & operator[](unsigned int); + * const value_type & operator[](unsigned int) const; + * }; + * @endcode + * + * @relates ReorderedIndexView + */ + template + internal::ReorderedIndexView + reordered_index_view(T &t) + { +#ifdef DEAL_II_WITH_CXX11 + static_assert(0 <= index && index < rank, + "The specified index must lie within the range [0,rank)"); +#endif + + return internal::ReorderedIndexView(t); + } + + + namespace internal + { + // ------------------------------------------------------------------------- + // Forward declarations and type traits + // ------------------------------------------------------------------------- + + template class StoreIndex; + template class Identity; + template class Contract2; + + /** + * An internally used type trait to allow nested application of the + * function reordered_index_view(T &t). + * + * The problem is that when working with the actual tensorial types, we + * have to return subtensors by reference - but sometimes, especially + * for StoreIndex and ReorderedIndexView that return rvalues, we have + * to return by value. + */ + template + struct ReferenceType + { + typedef T &type; + }; + + template + struct ReferenceType > + { + typedef StoreIndex type; + }; + + template + struct ReferenceType > + { + typedef ReorderedIndexView type; + }; + + + /** + * An internally used type trait that strips StoreIndex<0, S> down to + * its actual return value. This is needed to end the recursion in + * StoreIndex as well as, to return something meaningful in case of a + * nested application of the index reordering classes. + */ + template + struct StripStoreIndex + { + typedef T type; + }; + + template + struct StripStoreIndex > + { + typedef typename StripStoreIndex::type type; + }; + + + // TODO: Is there a possibility ot just have the following block of + // explanation on an internal page in doxygen? If, yes. Doxygen + // wizards, your call! + + // ------------------------------------------------------------------------- + // Implemenation of helper classes for reordered_index_view + // ------------------------------------------------------------------------- + + // OK. This is utterly brutal template magic. Therefore, we will not + // comment on the individual internal helper classes, because this is + // of not much value, but explain the general recursion procedure. + // + // (In order of appearance) + // + // Our task is to reorder access to a tensor object where a specified + // index is moved to the end. Thus we want to construct an object + // reorderd out of a tensor where the + // following access patterns are equivalent: + // @code + // tensor [i_0]...[i_index-1][i_index][i_index+1]...[i_n] + // reordered [i_0]...[i_index_1][i_index+1]...[i_n][i_index] + // @endcode + // + // The first task is to get rid of the application of + // [i_0]...[i_index-1]. This is a classical recursion pattern - relay + // the task from to by accessing the + // subtensor object: + + template + class ReorderedIndexView + { + public: + ReorderedIndexView(typename ReferenceType::type t) : t_(t) {} + + typedef ReorderedIndexView::value_type> + value_type; + + // Recurse by applying index j directly: + inline + value_type operator[](unsigned int j) const + { + return value_type(t_[j]); + } + + private: + typename ReferenceType::type t_; + }; + + // At some point we hit the condition index == 0, i.e., the first index + // should be reordered to the end. + // + // At this point we cannot be lazy any more and have to start storing + // indices because we get them in the wrong order. The user supplies + // [i_0][i_1]...[i_{rank - 1}] + // but we have to call the subtensor object with + // [i_{rank - 1}[i_0][i_1]...[i_{rank-2}] + // + // So give up and relay the task to the StoreIndex class: + + template + class ReorderedIndexView<0, rank, T> + { + public: + ReorderedIndexView(typename ReferenceType::type t) : t_(t) {} + + typedef internal::StoreIndex > value_type; + + inline + value_type operator[](unsigned int j) const + { + return value_type(internal::Identity(t_), j); + } + + private: + typename ReferenceType::type t_; + }; + + // Here, Identity is a helper class to ground the recursion in + // StoreIndex. Its implementation is easy - we haven't stored any + // indices yet. So, we just provide a function apply that returns the + // application of an index j to the stored tensor t_: + + template + class Identity + { + public: + Identity(typename ReferenceType::type t) : t_(t) {} + + typedef typename ValueType::value_type return_type; + + inline + typename ReferenceType::type apply(unsigned int j) const + { + return t_[j]; + } + + private: + typename ReferenceType::type t_; + }; + + // StoreIndex is a class that stores an index recursively with every + // invocation of operator[](unsigned int j): We do this by recursively + // creating a new StoreIndex class of lower rank that stores the + // supplied index j and holds a copy of the current class (with all + // other stored indices). Again, we provide an apply member function + // that knows how to apply an index on the highest rank and all + // subsequently stored indices: + + template + class StoreIndex : private S + { + public: + StoreIndex(S s, int i) : S(s), i_(i) {} + + typedef StoreIndex > value_type; + + inline + value_type operator[](unsigned int j) const + { + return value_type(*this, j); + } + + typedef typename ValueType::value_type return_type; + + inline + typename ReferenceType::type apply(unsigned int j) const + { + return S::apply(j)[i_]; + } + + private: + const int i_; + }; + + // We can store indices until we hit rank == 0. Then, we have all + // necessary indices and it is time to ground the recursion. For this, + // StoreIndex is specialized for rank == 0. The specialization contains + // to conversion operators to reference type to access the underlying + // object. Just call apply(i_) on the StoreIndex object of rank 1: + template + class StoreIndex<0, S> : private S + { + public: + StoreIndex(S s, int i) : S(s), i_(i) {} + + // Strip nested StoreIndex<0, S2> objects and cast to tensor's value_type + typedef typename StripStoreIndex::type value_type; + + inline operator value_type &() const + { + return S::apply(i_); + } + + private: + const int i_; + }; + + // ------------------------------------------------------------------------- + + } /* namespace internal */ } /* namespace TensorAccessors */ DEAL_II_NAMESPACE_CLOSE -- 2.39.5