diff --git a/dune/matrix-vector/genericvectortools.hh b/dune/matrix-vector/genericvectortools.hh index aee7664e9c72e34ffe27301f9c65305844322f78..9b29545930305792d4b9c9bb242e284662c8d5fb 100644 --- a/dune/matrix-vector/genericvectortools.hh +++ b/dune/matrix-vector/genericvectortools.hh @@ -47,97 +47,12 @@ void truncate(Vector& v, const BitVector& tr) { Impl::ScalarSwitch<Vector>::truncate(v, tr); } -//! Resize vector recursively to size of given vector/matrix -// EXPERIMENTAL: Its implementation is still subject to change. -// 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 -// TODO these are not satifsying, because e.g. for custom matrices that -// do not specialiaze the MatrixTraits properly, the error messages are -// not very helpful. In particular they may trigger a misleading static_assert. -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, typename std::enable_if_t<not IsNumber<Vector>::value>> { - using This = ScalarSwitch<Vector>; static void writeBinary(std::ostream& s, const Vector& v) { Hybrid::forEach(v, [&s](auto&& vi) { Generic::writeBinary(s, vi); }); @@ -153,75 +68,13 @@ struct ScalarSwitch<Vector, 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) { - // TODO this assert may also be triggered when isMatrix/isVector do not - // work properly, e.g. due to a missing MatrixTraits specialization. - 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. - /// We use AlwaysTrue to silence errors about enable_if not providing a type when it does not depend - /// on the function template parameter. - // ... from a sparse/dynamic matrix - template <class SizeProviderRow, typename = std::enable_if_t<isMatrix<Vector>() and not isTupleOrDerived<Vector>() and AlwaysTrue<SizeProviderRow>::value>> - 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>() and AlwaysTrue<SizeProviderRow>::value>> - 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>() and AlwaysTrue<SizeProviderRow>::value>> - 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) { - // TODO add static assert that fails or warns if subSizeProvider is not a matrix for a matrix or a vector for a vector - return This::subSizeProvider(std::forward<SizeProviderRow>(row), PriorityTag<42>()); - } }; //! recursion anchors for scalars template <class Scalar> struct ScalarSwitch<Scalar, typename std::enable_if_t<IsNumber<Scalar>::value>> { + static void writeBinary(std::ostream& s, const Scalar& v) { s.write(reinterpret_cast<const char*>(&v), sizeof(Scalar)); } @@ -237,8 +90,6 @@ struct ScalarSwitch<Scalar, v = 0; } - template <class Resizeable> - static void resize(Resizeable&, const Scalar&) {} }; } // end namespace Impl