Add resize for generic vector resizing.

dune/matrix-vector/genericvectortools.hh

@@ -5,11 +5,18 @@
  *  \brief Various tools for working with istl vectors of arbitrary nesting depth
  */
 
+#include <bitset>
 #include <iostream>
+#include <vector>
+
+#include <dune/common/classname.hh>
+#include <dune/common/concept.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/hybridutilities.hh>
+#include <dune/common/typetraits.hh>
 
 #include <dune/matrix-vector/algorithm.hh>
+#include <dune/matrix-vector/matrixtraits.hh>
 
 //! \brief Various tools for working with istl vectors of arbitrary nesting depth
 namespace Dune {
@@ -40,8 +47,89 @@ void truncate(Vector& v, const BitVector& tr) {
   Impl::ScalarSwitch<Vector>::truncate(v, tr);
 }
 
+//! Resize vector recursively to size of given vector/matrix
+// Note: We need the rvalue reference here, because we might obtain a
+// temporary recursively, e.g. in std::bitset the operator[] yields a
+// std::bitset::reference temporary due to its specific implementation.
+template <class Vector, class SizeProvider>
+void resize(Vector&& v, const SizeProvider& sizeProvider) {
+  Impl::ScalarSwitch<SizeProvider>::resize(v, sizeProvider);
+}
+// Overload for setting from a std::bitset. We need this because it
+// does not provide the range access the current implementation relies on.
+template <class Vector, long unsigned int n>
+void resize(Vector& v, const std::bitset<n>&) {
+  FieldVector<double, n> sizeProvider;
+  Impl::ScalarSwitch<decltype(sizeProvider)>::resize(v, sizeProvider);
+}
+
 namespace Impl {
 
+/// helper concepts
+namespace Concept {
+struct HasResize {
+  template <class C>
+  auto require(C&& c) -> decltype(c.resize(0));
+};
+
+struct HasSize {
+  // Note: Using Hybrid::size instead could be more general here. Make
+  //       sure to update the getSize correspondingly if necessary.
+  template <class C>
+  auto require(C&& c) -> decltype(c.size());
+};
+
+struct HasN {
+  template <class C>
+  auto require(C&& c) -> decltype(c.N());
+};
+
+} // namespace Concept
+
+//! Extract the size from the SizeProvider with these getSize overloads.
+/// Note: We use PriorityTags to enforce trial order. Call getSize without
+/// PriorityTag via the priority tag forward helper for standard behaviour.
+// look for size method
+template <class SizeProvider,
+          typename = std::enable_if_t<models<Concept::HasSize, SizeProvider>()>>
+auto getSize(const SizeProvider& sizeProvider, PriorityTag<20>) {
+  return sizeProvider.size();
+}
+// look for N method (matrices should have this)
+template <class SizeProvider,
+          typename = std::enable_if_t<models<Concept::HasN, SizeProvider>()>>
+auto getSize(const SizeProvider& sizeProvider, PriorityTag<10>) {
+  return sizeProvider.N();
+}
+// fail if none of the above succeeds
+template <class SizeProvider>
+auto getSize(const SizeProvider&, PriorityTag<0>) {
+  DUNE_THROW(NotImplemented, "Size extraction not implemented.");
+  return -1;
+}
+// priority tag forward helper
+template <class SizeProvider>
+auto getSize(const SizeProvider& sizeProvider) {
+  return getSize(sizeProvider, PriorityTag<42>());
+}
+
+/// compile-time helper traits
+template <class M>
+constexpr auto isMatrix() {
+  return Std::bool_constant<MatrixTraits<M>::isMatrix>();
+}
+
+template <class V>
+constexpr auto isVector() {
+  // Note: At the time this comment is written, the only use of
+  // isVector is in an overload after checking for isMatrix with
+  // lower priority tag, so isNotMatrix should not be needed here.
+  // Nevertheless, we keep it for clarity.
+  constexpr bool isNotMatrix = not decltype(isMatrix<V>())::value;
+  constexpr bool hasSize = models<Concept::HasSize, V>();
+  return Std::bool_constant<hasSize and isNotMatrix>();
+}
+
 //! recursion helper for scalars nested in iterables (vectors)
 template <class Vector>
 struct ScalarSwitch<Vector,
@@ -61,6 +149,65 @@ struct ScalarSwitch<Vector,
     sparseRangeFor(
         v, [&tr](auto&& vi, auto&& i) { Generic::truncate(vi, tr[i]); });
   }
+
+  // Note: The scalarSwitch is with respect to the sizeProvider to avoid
+  //       failures due to differing nesting depth. It might even be
+  //       intentional to have a sizeProvider with less nesting than the
+  //       vector to be resized.
+  template <class Resizeable>
+  static void resize(Resizeable& v, const Vector& sizeProvider) {
+    static_assert(not IsNumber<Resizeable>::value, "SizeProvider seems to have nesting depth that is not applicable to given vector type.");
+
+    using namespace Hybrid;
+    auto s = getSize(sizeProvider);
+    ifElse(models<Concept::HasResize, Resizeable>(),
+           [&](auto&& id) { id(v).resize(s); },
+           [&](auto) {
+             if (size(v) != s) {
+               std::stringstream msg;
+               msg << "No resize method found when trying to resize " << className(v) << " with size " << size(v) 
+                   << " from " << className(sizeProvider) << " with size " << s << ". "
+                   << "Note that you cannot resize statically sized vectors." << std::endl;
+               DUNE_THROW(RangeError, msg.str());
+             }
+           });
+
+    // recursively resize (in this branch of ScalarSwitch, all size providers are nested)
+    sparseRangeFor(sizeProvider, [&](auto&& row, auto&& i) {
+      Generic::resize(v[i], This::subSizeProvider(row));
+    });
+  }
+
+protected:
+  //! Obtain the subSizeProvider from a sizeProvider row element with these overloads.
+  /// Note: We use PriorityTags to enforce trial order w.r.t. failure fallback.
+  /// Additionally we use it to circumvent ambiguous overloading.
+  /// Call subSizeProvider without PriorityTag via the priority tag forward helper for standard behaviour.
+  // ... from a sparse/dynamic matrix
+  template <class SizeProviderRow, typename = std::enable_if_t<isMatrix<Vector>() and not isTupleOrDerived<Vector>()>>
+  static decltype(auto) subSizeProvider(SizeProviderRow&& row, PriorityTag<21>) {
+    return *row.begin();
+  }
+  // .. from a tuple like/static matrix
+  template <class SizeProviderRow, typename = std::enable_if_t<isMatrix<Vector>() and isTupleOrDerived<Vector>()>>
+  static decltype(auto) subSizeProvider(SizeProviderRow&& row, PriorityTag<20>) {
+    return row[Indices::_0];
+  }
+  // ... from a vector
+  template <class SizeProviderRow, typename = std::enable_if_t<isVector<Vector>()>>
+  static decltype(auto) subSizeProvider(SizeProviderRow&& row, PriorityTag<10>) {
+    return row;
+  }
+  // fail when none of the above can be applied to avoid unexpected resize behaviour
+  template <class SizeProviderRow>
+  static void subSizeProvider(SizeProviderRow&&, PriorityTag<0>) {
+    static_assert(AlwaysFalse<Vector>::value, "MatrixVector::Generic::resize does not seem to work properly on given SizeProvider type. If it is a matrix, maybe you need to explicitly inject Dune::MatrixVector::MatrixTraits<SizeProvider>::isMatrix (see dune-matrix-vector/dune/matrix-vector/matrixtraits.hh).");
+  }
+  // priority tag forward helper
+  template <class SizeProviderRow>
+  static decltype(auto) subSizeProvider(SizeProviderRow&& row) {
+    return This::subSizeProvider(std::forward<SizeProviderRow>(row), PriorityTag<42>());
+  }
 };
 
 //! recursion anchors for scalars
@@ -81,6 +228,9 @@ struct ScalarSwitch<Scalar,
     if (tr)
       v = 0;
   }
+
+  template <class Resizeable>
+  static void resize(Resizeable&, const Scalar&) {}
 };
 
 } // end namespace Impl