diff --git a/dune/matrix-vector/resize.hh b/dune/matrix-vector/resize.hh new file mode 100644 index 0000000000000000000000000000000000000000..41ec988465709636321835f112575b7333cf9e94 --- /dev/null +++ b/dune/matrix-vector/resize.hh @@ -0,0 +1,130 @@ +#ifndef DUNE_MATRIX_VECTOR_RESIZE_HH +#define DUNE_MATRIX_VECTOR_RESIZE_HH + +#include <dune/common/fmatrix.hh> +#include <dune/common/fvector.hh> +#include <dune/common/hybridutilities.hh> +#include <dune/common/typetraits.hh> +#include <dune/common/typeutilities.hh> + +#include <dune/istl/bvector.hh> +#include <dune/istl/bcrsmatrix.hh> + +#include <dune/matrix-vector/algorithm.hh> +#include <dune/matrix-vector/concepts.hh> +#include <dune/matrix-vector/matrixtraits.hh> +#include <dune/matrix-vector/traitutilities.hh> + +namespace Dune { +namespace MatrixVector { + +// forward declaration +template <class Vector, class SizeProvider> +void resize(Vector& v, const SizeProvider& sizeProvider); + +namespace Impl { + +//! checks for resizeability and resizes +template <class Vector> +void tryResize(Vector& v, size_t size) { + Hybrid::ifElse( + models<Concept::HasResize, Vector>(), + [&](auto&& id) { id(v).resize(size); }, + [&](auto) { + if (Hybrid::size(v) != size) { + std::stringstream msg; + msg << "No resize method found when trying to resize " << className(v) + << " with size " << Hybrid::size(v) << "to size " << size << ". " + << "Note that you cannot resize statically sized vectors." + << std::endl; + DUNE_THROW(RangeError, msg.str()); + } + }); +} + +/// Implementations for recursive resize. +/// Note: We incorporate priority tags to resolve special cases in well defined order +/// (e.g. resizing a std::bitset from another std::bitset), but also to resolve SFINAE +/// overload ambiguity due to the Enable* terms. +// Treat a size providing std::bitset as if it was a FieldVector of same size. +template <class Vector, unsigned long int n> +void resize(Vector& v, const std::bitset<n>&, PriorityTag<11>) { + Dune::MatrixVector::resize(v, FieldVector<double, n>()); +} +// Resizing an std::bitset causes is a failure if static and desired size do not coincide. +template <unsigned long int n, class SizeProvider> +void resize(std::bitset<n>& v, const SizeProvider& sizeProvider, PriorityTag<10>) { + using namespace Dune::Hybrid; + ifElse(models<Concept::HasSize, SizeProvider>(), [&] (auto) { + if(sizeProvider.size() != n) + DUNE_THROW(RangeError, "Cannot resize statically sized std::bitset<n> to different value."); + }); + ifElse(not models<Concept::HasSize, SizeProvider>() and + models<Concept::HasN, SizeProvider>(), [&] (auto) { + if(sizeProvider.size() != n) + DUNE_THROW(RangeError, "Cannot resize statically sized std::bitset<n> to different value."); + }); + ifElse(not models<Concept::HasSize, SizeProvider()>() and + not models<Concept::HasN, SizeProvider>(), [&] (auto id) { + // Note: This could be a static_assert but I dont't know how to make it compile. + DUNE_THROW(NotImplemented, "Size provider does neither provide a method size() nor N()."); + }); +} + +// generic resize recursion anchor for scalars +template <class Vector, class SizeProvider, typename = EnableScalar<SizeProvider>> +void resize(Vector&, const SizeProvider&, PriorityTag<3>) { + static_assert(isScalar<Vector>(), + "Vector nesting depth exceeds size provider nesting."); +} +// generic resize from vector +template <class Vector, class SizeProvider, typename = EnableVector<SizeProvider>> +void resize(Vector& v, const SizeProvider& sizeProvider, PriorityTag<2>) { + using namespace Hybrid; + // resize + static_assert(models<Concept::HasSize, SizeProvider>(), + "Size provider was considered a vetor but does not provide a method size()."); + tryResize(v, sizeProvider.size()); + // resize recursively + static_assert(isSparseRangeIterable<SizeProvider>(), + "Cannot resize recursively due to missing implementation for size provider type."); + sparseRangeFor(sizeProvider, [&](auto&& subSizeProvider, auto&& i) { + Dune::MatrixVector::resize(v[i], subSizeProvider); + }); +} +// generic resize from matrix +template <class Vector, class SizeProvider, typename = EnableMatrix<SizeProvider>> +void resize(Vector& v, const SizeProvider& sizeProvider, PriorityTag<1>) { + using namespace Dune::Hybrid; + // resize + static_assert(models<Concept::HasN, SizeProvider>(), + "Size provider was considered a matrix but does not provide a method N()."); + tryResize(v, sizeProvider.N()); + // recursive resize + static_assert(isSparseRangeIterable<SizeProvider>(), + "Cannot resize recursively due to missing implementation for size provider type."); + sparseRangeFor(sizeProvider, [&](auto&& row, auto&& i) { + static_assert(isSparseRangeIterable<decltype(row)>(), + "Cannot resize recursively due to missing implementation for size provider row type."); + sparseRangeFirst(row, [&](auto&& subSizeProvider) { + Dune::MatrixVector::resize(v[i], subSizeProvider); + }); + }); +} +// fallback compile-time error +template <class Vector, class SizeProvider> +void resize(Vector&, const SizeProvider&, PriorityTag<0>) { + static_assert(AlwaysFalse<Vector>::value, "No resize method found."); +} + +} // end namespace Impl + +template <class Vector, class SizeProvider> +void resize(Vector& v, const SizeProvider& sizeProvider) { + Impl::resize(v, sizeProvider, PriorityTag<42>()); +} + +} // end namespace MatrixVector +} // end namespace Dune + +#endif // DUNE_MATRIX_VECTOR_RESIZE_HH