]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement TensorAccessors::reordered_index_view
authorMatthias Maier <tamiko@43-1.org>
Fri, 11 Sep 2015 19:41:54 +0000 (14:41 -0500)
committerMatthias Maier <tamiko@43-1.org>
Fri, 11 Sep 2015 19:41:54 +0000 (14:41 -0500)
include/deal.II/base/tensor_accessors.h

index 22e9fcf6f338cdd5f54934ace18ab7f1dba3a1b1..9d4e405a7fed2244f43ea8b0a467ee19379e4c53 100644 (file)
@@ -68,8 +68,353 @@ DEAL_II_NAMESPACE_OPEN
  */
 namespace TensorAccessors
 {
+  // forward declarations
+  namespace internal
+  {
+    template <int index, int rank, typename T> 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
+   *  - <code>T::value_type</code> if T is a tensorial class providing a
+   *    typedef <code>value_type</code> and does not have a const qualifier.
+   *  - <code>const T::value_type</code> if T is a tensorial class
+   *    providing a typedef <code>value_type</code> and does have a const
+   *    qualifier.
+   *  - <code>const T::value_type</code> if T is a tensorial class
+   *    providing a typedef <code>value_type</code> and does have a const
+   *    qualifier.
+   *  - <code>A</code> if T is of array type <code>A[...]</code>
+   *  - <code>const A</code> if T is of array type <code>A[...]</code> and
+   *    does have a const qualifier.
+   */
+  template <typename T>
+  struct ValueType
+  {
+    typedef typename T::value_type value_type;
+  };
+
+  template <typename T>
+  struct ValueType<const T>
+  {
+    typedef const typename T::value_type value_type;
+  };
+
+  template <typename T, std::size_t N>
+  struct ValueType<T[N]>
+  {
+    typedef T value_type;
+  };
+
+  template <typename T, std::size_t N>
+  struct ValueType<const T[N]>
+  {
+    typedef const T value_type;
+  };
+
+
+  /**
+   * This class provides a local typedef @p value_type that is equal to
+   * the typedef <code>value_type</code> 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<rank, T>::value_type
+   * will also be const.
+   */
+  template <int deref_steps, typename T>
+  struct ReturnType
+  {
+    typedef typename ReturnType<deref_steps - 1, typename ValueType<T>::value_type>::value_type value_type;
+  };
+
+  template <typename T>
+  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
+   * <code>operator[]()</code> 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 <code>operator[](unsigned int)</code> and a
+   * typedef <code>value_type</code> 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 <code>value_type</code> and an index operator
+   * <code>operator[]()</code> that returns a (const or non-const)
+   * reference of <code>value_type</code>:
+   * @code
+   *   class T
+   *   {
+   *     typedef ... value_type
+   *     value_type & operator[](unsigned int);
+   *     const value_type & operator[](unsigned int) const;
+   *   };
+   * @endcode
+   *
+   * @relates ReorderedIndexView
+   */
+  template <int index, int rank, typename T>
+  internal::ReorderedIndexView<index, rank, T>
+  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<index, rank, T>(t);
+  }
+
+
+  namespace internal
+  {
+    // -------------------------------------------------------------------------
+    // Forward declarations and type traits
+    // -------------------------------------------------------------------------
+
+    template <int rank, typename S> class StoreIndex;
+    template <typename T> class Identity;
+    template <int no_contr, int dim> 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<typename T>
+    struct ReferenceType
+    {
+      typedef T &type;
+    };
+
+    template <int rank, typename S>
+    struct ReferenceType<StoreIndex<rank, S> >
+    {
+      typedef StoreIndex<rank, S> type;
+    };
+
+    template <int index, int rank, typename T>
+    struct ReferenceType<ReorderedIndexView<index, rank, T> >
+    {
+      typedef ReorderedIndexView<index, rank, T> 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<typename T>
+    struct StripStoreIndex
+    {
+      typedef T type;
+    };
+
+    template <typename S>
+    struct StripStoreIndex<StoreIndex<0, S> >
+    {
+      typedef typename StripStoreIndex<typename S::return_type>::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
+    // <code>reorderd</code> out of a <code>tensor</code> 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 <index, rank> to <index-1, rank-1> by accessing the
+    // subtensor object:
+
+    template <int index, int rank, typename T>
+    class ReorderedIndexView
+    {
+    public:
+      ReorderedIndexView(typename ReferenceType<T>::type t) : t_(t) {}
+
+      typedef ReorderedIndexView<index - 1, rank - 1, typename ValueType<T>::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<T>::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 <int rank, typename T>
+    class ReorderedIndexView<0, rank, T>
+    {
+    public:
+      ReorderedIndexView(typename ReferenceType<T>::type t) : t_(t) {}
+
+      typedef internal::StoreIndex<rank - 1, internal::Identity<T> > value_type;
+
+      inline
+      value_type operator[](unsigned int j) const
+      {
+        return value_type(internal::Identity<T>(t_), j);
+      }
+
+    private:
+      typename ReferenceType<T>::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 <typename T>
+    class Identity
+    {
+    public:
+      Identity(typename ReferenceType<T>::type t) : t_(t) {}
+
+      typedef typename ValueType<T>::value_type return_type;
+
+      inline
+      typename ReferenceType<return_type>::type apply(unsigned int j) const
+      {
+        return t_[j];
+      }
+
+    private:
+      typename ReferenceType<T>::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 <int rank, typename S>
+    class StoreIndex : private S
+    {
+    public:
+      StoreIndex(S s, int i) : S(s), i_(i) {}
+
+      typedef StoreIndex<rank - 1, StoreIndex<rank, S> > value_type;
+
+      inline
+      value_type operator[](unsigned int j) const
+      {
+        return value_type(*this, j);
+      }
+
+      typedef typename ValueType<typename S::return_type>::value_type return_type;
+
+      inline
+      typename ReferenceType<return_type>::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 <typename S>
+    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<typename S::return_type>::type value_type;
+
+      inline operator value_type &() const
+      {
+        return S::apply(i_);
+      }
+
+    private:
+      const int i_;
+    };
+
+    // -------------------------------------------------------------------------
+
+  } /* namespace internal */
 } /* namespace TensorAccessors */
 
 DEAL_II_NAMESPACE_CLOSE

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.