Bugfix: Fix return value of StoredIndex<1, S>

The problem is subtle. TensorAccessors::reordered_index_view returns an
internal object that behaves like a tensor.

The problem is that fully accessing the tensorial object, i.e., applying
operator[](unsigned int) rank times returned an object of type

  StoreIndex<0, [...]>

instead of type Number of the original tensor. The problem is that the
compiler has no way of knowing whether we want to cast this object to
Number during type matching, thus something like

  std::min( StoreIndex<0, [...]>(...), (double)4.);

fails.

diff --git a/include/deal.II/base/tensor_accessors.h b/include/deal.II/base/tensor_accessors.h
index 97a9478553f3dc995fb4b4b705fbda83dfce2c39..a80057dd3760a5c43d41be6f0ff8dd2b89c92ab4 100644 (file)
--- a/include/deal.II/base/tensor_accessors.h
+++ b/include/deal.II/base/tensor_accessors.h
@@ -355,25 +355,6 @@ namespace TensorAccessors
     };
 
 
-    /**
-     * 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!
@@ -422,8 +403,8 @@ namespace TensorAccessors
       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 some point we hit the condition index == 0 and rank > 1, 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
@@ -439,12 +420,32 @@ namespace TensorAccessors
     public:
       ReorderedIndexView(typename ReferenceType<T>::type t) : t_(t) {}
 
-      typedef internal::StoreIndex<rank - 1, internal::Identity<T> > value_type;
+      typedef StoreIndex<rank - 1, internal::Identity<T> > value_type;
 
       inline
       value_type operator[](unsigned int j) const
       {
-        return value_type(internal::Identity<T>(t_), j);
+        return value_type(Identity<T>(t_), j);
+      }
+
+    private:
+      typename ReferenceType<T>::type t_;
+    };
+
+    // Sometimes, we're lucky and don't have to do anything. In this case
+    // just return the original tensor.
+
+    template <typename T>
+    class ReorderedIndexView<0, 1, T>
+    {
+    public:
+      ReorderedIndexView(typename ReferenceType<T>::type t) : t_(t) {}
+
+      typedef typename ReferenceType<typename ValueType<T>::value_type>::type value_type;
+
+      inline value_type operator[](unsigned int j) const
+      {
+        return t_[j];
       }
 
     private:
@@ -483,7 +484,7 @@ namespace TensorAccessors
     // subsequently stored indices:
 
     template <int rank, typename S>
-    class StoreIndex : private S
+    class StoreIndex : public S
     {
     public:
       StoreIndex(S s, int i) : S(s), i_(i) {}
@@ -508,23 +509,22 @@ namespace TensorAccessors
       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:
+    // We have to store indices until we hit rank == 1. Then, upon the next
+    // invocation of operator[](unsigned int j) we have all necessary
+    // information available to return the actual object.
+
     template <typename S>
-    class StoreIndex<0, S> : private S
+    class StoreIndex<1, S> : public 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;
+      typedef typename ValueType<typename S::return_type>::value_type return_type;
+      typedef return_type value_type;
 
-      inline operator value_type &() const
+      inline return_type &operator[](unsigned int j) const
       {
-        return S::apply(i_);
+        return S::apply(j)[i_];
       }
 
     private:

diff --git a/tests/base/tensor_accessors_01.cc b/tests/base/tensor_accessors_01.cc
index 8e51c76dd6bbe73a2c2360d44f2447540143e0ae..10342da329e7668a6f71ee2837b6af614f84b7e4 100644 (file)
--- a/tests/base/tensor_accessors_01.cc
+++ b/tests/base/tensor_accessors_01.cc
@@ -127,4 +127,13 @@ int main()
     foo2 = TensorAccessors::reordered_index_view<2, 5>(t_ref);
     deallog << foo2[0][1][0][1][2] << std::endl;
   }
+
+  {
+    // Is it possible to call the simplest case (where we have to do
+    // absolutely nothing?
+
+    Tensor<1, 3, int> t;
+    TensorAccessors::internal::ReorderedIndexView<0, 1, Tensor<1, 3, int> > // auto ...
+    foo = TensorAccessors::reordered_index_view<0, 1>(t);
+  }
 }