diff --git a/dune/matrix-vector/genericvectortools.hh b/dune/matrix-vector/genericvectortools.hh index 27aa99594889e8276b2154a76fc05063ae877c63..9c15e8a4aa5664085ecf37c5180da3f42716334a 100644 --- a/dune/matrix-vector/genericvectortools.hh +++ b/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