Implement TensorAccessors::extract

diff --git a/include/deal.II/base/tensor_accessors.h b/include/deal.II/base/tensor_accessors.h
index 9d4e405a7fed2244f43ea8b0a467ee19379e4c53..80cd779e2dc61f9830d8b1cfd38e2b75c2557f20 100644 (file)
--- a/include/deal.II/base/tensor_accessors.h
+++ b/include/deal.II/base/tensor_accessors.h
@@ -72,6 +72,7 @@ namespace TensorAccessors
   namespace internal
   {
     template <int index, int rank, typename T> class ReorderedIndexView;
+    template <int position, int rank> struct ExtractHelper;
   }
 
 
@@ -198,6 +199,34 @@ namespace TensorAccessors
   }
 
 
+  /**
+   * Return a reference (const or non-const) to a subobject of a tensorial
+   * object @p t of type @p T, as described by an array type @p ArrayType
+   * object @p indices. For example: @code
+   *   Tensor<5, dim> tensor;
+   *   TableIndices<5> indices (0, 1, 2, 3, 4);
+   *   TensorAccessors::extract(tensor, indices) = 42;
+   * @endcode
+   * This is equivalent to <code>tensor[0][1][2][3][4] = 42.</code>.
+   *
+   * @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>. Further, its tensorial rank must
+   * be equal or greater than @p rank
+   *
+   * @tparam ArrayType An array like object, such as std::array, or
+   * dealii::TableIndices  that stores at least @p rank indices that can be
+   * accessed via operator[]().
+   */
+  template<int rank, typename T, typename ArrayType> typename
+  ReturnType<rank, T>::value_type &
+  extract(T &t, const ArrayType &indices)
+  {
+    return internal::ExtractHelper<0, rank>::template extract<T, ArrayType>(t, indices);
+  }
+
+
   namespace internal
   {
     // -------------------------------------------------------------------------
@@ -412,6 +441,41 @@ namespace TensorAccessors
       const int i_;
     };
 
+
+    // -------------------------------------------------------------------------
+    // Implemenation of helper classes for extract
+    // -------------------------------------------------------------------------
+
+    // Straightforward recursion implemented by specializing ExtractHelper
+    // for position == rank. We use the type trait ReturnType<rank, T> to
+    // have an idea what the final type will be.
+    template<int position, int rank>
+    struct ExtractHelper
+    {
+      template<typename T, typename ArrayType>
+      inline static
+      typename ReturnType<rank - position, T>::value_type &
+      extract(T &t, const ArrayType &indices)
+      {
+        return ExtractHelper<position + 1, rank>::
+               template extract<typename ValueType<T>::value_type, ArrayType>
+        (t[indices[position]], indices);
+      }
+    };
+
+    // For dimension == rank there is nothing to extract, just return the
+    // object.
+    template<int rank>
+    struct ExtractHelper<rank, rank>
+    {
+      template<typename T, typename ArrayType>
+      inline static
+      T &extract(T &t, const ArrayType &indices)
+      {
+        return t;
+      }
+    };
+
     // -------------------------------------------------------------------------
 
   } /* namespace internal */