Loading...
Searching...
No Matches
istlvectorbackend.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_ISTLVECTORBACKEND_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_ISTLVECTORBACKEND_HH
#include <cstddef>
#include <utility>
#include <type_traits>
#include <dune/common/std/type_traits.hh>
#include <dune/common/indices.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/concept.hh>
#include <dune/functions/common/indexaccess.hh>
#include <dune/functions/functionspacebases/concepts.hh>
namespace Dune {
namespace Functions {
namespace Impl {
template<class V,
std::enable_if_t<not Dune::models<Imp::Concept::HasStaticIndexAccess, V>() , int> = 0>
auto fieldTypes(V&& /*v*/)
{
return TypeList<V>{};
}
template<class V,
std::enable_if_t<Dune::models<Imp::Concept::HasStaticIndexAccess, V>(), int> = 0>
auto fieldTypes(V&& v)
{
if constexpr (Dune::models<Imp::Concept::HasDynamicIndexAccess<std::size_t>, V>())
return fieldTypes(v[std::size_t{0}]);
else
{
auto indexRange = typename decltype(range(Hybrid::size(v)))::integer_sequence();
return unpackIntegerSequence([&](auto... i) {
return uniqueTypeList(std::tuple_cat(fieldTypes(v[i])...));
}, indexRange);
}
}
} // namespace Impl
template<class V>
constexpr auto fieldTypes()
{
return decltype(Impl::fieldTypes(std::declval<V>())){};
}
template<class V>
constexpr bool hasUniqueFieldType()
{
return std::tuple_size_v<std::decay_t<decltype(fieldTypes<V>())>> ==1;
}
namespace Impl {
/*
* \brief A wrapper providing multi-index access to vector entries
*
* The wrapped vector type should be an istl like random
* access container providing operator[] and size() methods.
* For classical containers this should support indices
* of type std::size_t. For multi-type containers indices
* of the form Dune::index_constant<i> should be supported
* while size() should be a static constexpr method.
*
* When resolving multi-indices the backend appends indices
* using operator[] as long as the result is not a scalar.
* If this exhausts the digits of the multi-index, additional
* zeros are appended.
*
* \tparam V Type of the raw wrapper vector
*/
template<class V>
class ISTLVectorBackend
{
// Template aliases for using detection idiom.
template<class C>
using dynamicIndexAccess_t = decltype(std::declval<C>()[0]);
template<class C>
using staticIndexAccess_t = decltype(std::declval<C>()[Dune::Indices::_0]);
template<class C>
using resizeMethod_t = decltype(std::declval<C>().resize(0));
// Short cuts for feature detection
template<class C>
using hasDynamicIndexAccess = Dune::Std::is_detected<dynamicIndexAccess_t, std::remove_reference_t<C>>;
template<class C>
using hasStaticIndexAccess = Dune::Std::is_detected<staticIndexAccess_t, std::remove_reference_t<C>>;
template<class C>
using hasResizeMethod = Dune::Std::is_detected<resizeMethod_t, std::remove_reference_t<C>>;
template<class C>
using isDynamicVector = Dune::Std::is_detected<dynamicIndexAccess_t, std::remove_reference_t<C>>;
template<class C>
using isStaticVector = std::bool_constant<
Dune::Std::is_detected_v<staticIndexAccess_t, std::remove_reference_t<C>>
and not Dune::Std::is_detected_v<dynamicIndexAccess_t, std::remove_reference_t<C>>>;
template<class C>
using isScalar = std::bool_constant<not Dune::Std::is_detected_v<staticIndexAccess_t, std::remove_reference_t<C>>>;
template<class C>
using isVector = std::bool_constant<Dune::Std::is_detected_v<staticIndexAccess_t, std::remove_reference_t<C>>>;
template<class... Args>
static void forwardToResize(Args&&... args)
{
resize(std::forward<Args>(args)...);
}
template<class C, class SizeProvider,
std::enable_if_t<hasResizeMethod<C>::value, int> = 0>
static void resize(C&& c, const SizeProvider& sizeProvider, typename SizeProvider::SizePrefix prefix)
{
auto size = sizeProvider.size(prefix);
if (size==0)
{
// If size==0 this prefix refers to a single coefficient c.
// But being in this overload means that c is not a scalar
// because is has a resize() method. Since operator[] deliberately
// supports implicit padding of multi-indices by as many
// [0]'s as needed to resolve a scalar entry, we should also
// except a non-scalar c here. However, this requires that
// we silently believe that whatever size c already has is
// intended by the user. The only exception is c.size()==0
// which is not acceptable but we also cannot know the desired size.
if (c.size()==0)
DUNE_THROW(RangeError, "The vector entry v[" << prefix << "] should refer to a "
<< "scalar coefficient, but is a dynamically sized vector of size==0");
else
// Accept non-zero sized coefficients to avoid that resize(basis)
// fails for a vector that works with operator[] and already
// has the appropriate size.
return;
}
c.resize(size);
prefix.push_back(0);
for(std::size_t i=0; i<size; ++i)
{
prefix.back() = i;
resize(c[i], sizeProvider, prefix);
}
}
template<class C, class SizeProvider,
std::enable_if_t<not hasResizeMethod<C>::value, int> = 0,
std::enable_if_t<isVector<C>::value, int> = 0>
static void resize(C&& c, const SizeProvider& sizeProvider, typename SizeProvider::SizePrefix prefix)
{
auto size = sizeProvider.size(prefix);
// If size == 0 there's nothing to do:
// We can't resize c and it's already
// large enough anyway.
if (size == 0)
return;
// If size>0 but c does not have the appropriate
// size we throw an exception.
//
// We could perhaps also allow c.size()>size.
// But then looping the loop below gets complicated:
// We're not allowed to loop until c.size(). But
// we also cannot use size for termination,
// because this fails if c is a static vector.
if (c.size() != size)
DUNE_THROW(RangeError, "Can't resize non-resizable entry v[" << prefix << "] of size " << c.size() << " to size(" << prefix << ")=" << size);
// Recursively resize all entries of c now.
using namespace Dune::Hybrid;
prefix.push_back(0);
forEach(integralRange(Hybrid::size(c)), [&](auto&& i) {
prefix.back() = i;
// Here we'd simply like to call resize(c[i], sizeProvider, prefix);
// but even gcc-7 does not except this bus reports
// "error: ‘this’ was not captured for this lambda function"
// although there's no 'this' because we're in a static method.
// Bypassing this by and additional method that does perfect
// forwarding allows to workaround this.
ISTLVectorBackend<V>::forwardToResize(c[i], sizeProvider, prefix);
});
}
template<class C, class SizeProvider,
std::enable_if_t<not hasResizeMethod<C>::value, int> = 0,
std::enable_if_t<isScalar<C>::value, int> = 0>
static void resize(C&&, const SizeProvider& sizeProvider, typename SizeProvider::SizePrefix prefix)
{
auto size = sizeProvider.size(prefix);
if (size != 0)
DUNE_THROW(RangeError, "Can't resize scalar vector entry v[" << prefix << "] to size(" << prefix << ")=" << size);
}
template<class C, class T,
std::enable_if_t<std::is_assignable_v<C&,T>, int> = 0>
void recursiveAssign(C& c, const T& t)
{
c = t;
}
template<class C, class T,
std::enable_if_t<not std::is_assignable_v<C&,T>, int> = 0>
void recursiveAssign(C& c, const T& t)
{
Dune::Hybrid::forEach(c, [&](auto&& ci) {
recursiveAssign(ci, t);
});
}
public:
using Vector = V;
ISTLVectorBackend(Vector& vector) :
vector_(&vector)
{}
template<class SizeProvider>
void resize(const SizeProvider& sizeProvider)
{
auto prefix = typename SizeProvider::SizePrefix();
prefix.resize(0);
resize(*vector_, sizeProvider, prefix);
}
template<class MultiIndex>
decltype(auto) operator[](const MultiIndex& index) const
{
return resolveDynamicMultiIndex(*vector_, index);
}
template<class MultiIndex>
decltype(auto) operator[](const MultiIndex& index)
{
return resolveDynamicMultiIndex(*vector_, index);
}
template<typename T>
void operator= (const T& other)
{
recursiveAssign(vector(), other);
}
template<typename T>
void operator= (const ISTLVectorBackend<T>& other)
{
vector() = other.vector();
}
const Vector& vector() const
{
return *vector_;
}
Vector& vector()
{
return *vector_;
}
private:
Vector* vector_;
};
} // end namespace Impl
template<class Vector>
auto istlVectorBackend(Vector& v)
{
static_assert(hasUniqueFieldType<Vector&>(), "Vector type passed to istlVectorBackend() does not have a unique field type.");
return Impl::ISTLVectorBackend<Vector>(v);
}
template<class Vector>
auto istlVectorBackend(const Vector& v)
{
static_assert(hasUniqueFieldType<const Vector&>(), "Vector type passed to istlVectorBackend() does not have a unique field type.");
return Impl::ISTLVectorBackend<const Vector>(v);
}
} // namespace Functions
} // namespace Dune
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_ISTLVECTORBACKEND_HH
indexaccess.hh
Dune::Functions::istlVectorBackend
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition istlvectorbackend.hh:350
Dune::Functions::resolveDynamicMultiIndex
constexpr decltype(auto) resolveDynamicMultiIndex(C &&c, const MultiIndex &multiIndex, const IsFinal &isFinal)
Provide multi-index access by chaining operator[].
Definition indexaccess.hh:377
Dune
Definition monomialset.hh:19
Dune::Functions
Definition monomialset.hh:19
Dune::Functions::fieldTypes
constexpr auto fieldTypes()
Generate list of field types in container.
Definition istlvectorbackend.hh:67
Dune::Functions::hasUniqueFieldType
constexpr bool hasUniqueFieldType()
Check if container has a unique field type.
Definition istlvectorbackend.hh:78
concepts.hh