DUNE PDELab (2.8)

istlvectorbackend.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_ISTLVECTORBACKEND_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_ISTLVECTORBACKEND_HH
5
6#include <cstddef>
7#include <utility>
8#include <type_traits>
9
10#include <dune/common/std/type_traits.hh>
11#include <dune/common/indices.hh>
12#include <dune/common/hybridutilities.hh>
14
15#include <dune/functions/common/indexaccess.hh>
16#include <dune/functions/functionspacebases/concepts.hh>
17
18
19namespace Dune {
20namespace Functions {
21
22namespace Impl {
23
24template<class V,
25 std::enable_if_t<not Dune::models<Imp::Concept::HasStaticIndexAccess, V>() , int> = 0>
26auto fieldTypes(V&& v)
27{
28 return TypeList<V>{};
29}
30
31template<class V,
32 std::enable_if_t<Dune::models<Imp::Concept::HasStaticIndexAccess, V>(), int> = 0>
33auto fieldTypes(V&& v)
34{
35 if constexpr (Dune::models<Imp::Concept::HasDynamicIndexAccess<std::size_t>, V>())
36 return fieldTypes(v[std::size_t{0}]);
37 else
38 {
39 auto indexRange = typename decltype(range(Hybrid::size(v)))::integer_sequence();
40 return unpackIntegerSequence([&](auto... i) {
41 return uniqueTypeList(std::tuple_cat(fieldTypes(v[i])...));
42 }, indexRange);
43 }
44}
45
46} // namespace Impl
47
48
49
62template<class V>
63constexpr auto fieldTypes()
64{
65 return decltype(Impl::fieldTypes(std::declval<V>())){};
66}
67
73template<class V>
74constexpr bool hasUniqueFieldType()
75{
76 return std::tuple_size<std::decay_t<decltype(fieldTypes<V>())>>::value==1;
77}
78
79
80
81namespace Impl {
82
83/*
84 * \brief A wrapper providing multi-index access to vector entries
85 *
86 * The wrapped vector type should be an istl like random
87 * access container providing operator[] and size() methods.
88 * For classical containers this should support indices
89 * of type std::size_t. For multi-type containers indices
90 * of the form Dune::index_constant<i> should be supported
91 * while size() should be a static constexpr method.
92 *
93 * When resolving multi-indices the backend appends indices
94 * using operator[] as long as the result is not a scalar.
95 * If this exhausts the digits of the multi-index, additional
96 * zero`s are appended.
97 *
98 * \tparam V Type of the raw wrapper vector
99 */
100template<class V>
101class ISTLVectorBackend
102{
103
104 // Template aliases for using detection idiom.
105 template<class C>
106 using dynamicIndexAccess_t = decltype(std::declval<C>()[0]);
107
108 template<class C>
109 using staticIndexAccess_t = decltype(std::declval<C>()[Dune::Indices::_0]);
110
111 template<class C>
112 using resizeMethod_t = decltype(std::declval<C>().resize(0));
113
114
115
116 // Short cuts for feature detection
117 template<class C>
119
120 template<class C>
122
123 template<class C>
125
126 template<class C>
128
129 template<class C>
130 using isStaticVector = Dune::Std::bool_constant<
131 Dune::Std::is_detected_v<staticIndexAccess_t, std::remove_reference_t<C>>
132 and not Dune::Std::is_detected_v<dynamicIndexAccess_t, std::remove_reference_t<C>>>;
133
134 template<class C>
135 using isScalar = Dune::Std::bool_constant<not Dune::Std::is_detected_v<staticIndexAccess_t, std::remove_reference_t<C>>>;
136
137 template<class C>
138 using isVector = Dune::Std::bool_constant<Dune::Std::is_detected_v<staticIndexAccess_t, std::remove_reference_t<C>>>;
139
140
141
142 template<class... Args>
143 static void forwardToResize(Args&&... args)
144 {
145 resize(std::forward<Args>(args)...);
146 }
147
148
149 template<class C, class SizeProvider,
150 std::enable_if_t<hasResizeMethod<C>::value, int> = 0>
151 static void resize(C&& c, const SizeProvider& sizeProvider, typename SizeProvider::SizePrefix prefix)
152 {
153 auto size = sizeProvider.size(prefix);
154 if (size==0)
155 {
156 // If size==0 this prefix refers to a single coefficient c.
157 // But being in this overload means that c is not a scalar
158 // because is has a resize() method. Since operator[] deliberately
159 // supports implicit padding of multi-indices by as many
160 // [0]'s as needed to resolve a scalar entry, we should also
161 // except a non-scalar c here. However, this requires that
162 // we silently believe that whatever size c already has is
163 // intended by the user. The only exception is c.size()==0
164 // which is not acceptable but we also cannot know the desired size.
165 if (c.size()==0)
166 DUNE_THROW(RangeError, "The vector entry v[" << prefix << "] should refer to a "
167 << "scalar coefficient, but is a dynamically sized vector of size==0");
168 else
169 // Accept non-zero sized coefficients to avoid that resize(basis)
170 // fails for a vector that works with operator[] and already
171 // has the appropriate size.
172 return;
173 }
174 c.resize(size);
175 prefix.push_back(0);
176 for(std::size_t i=0; i<size; ++i)
177 {
178 prefix.back() = i;
179 resize(c[i], sizeProvider, prefix);
180 }
181 }
182
183 template<class C, class SizeProvider,
184 std::enable_if_t<not hasResizeMethod<C>::value, int> = 0,
185 std::enable_if_t<isVector<C>::value, int> = 0>
186 static void resize(C&& c, const SizeProvider& sizeProvider, typename SizeProvider::SizePrefix prefix)
187 {
188 auto size = sizeProvider.size(prefix);
189 // If size == 0 there's nothing to do:
190 // We can't resize c and it's already
191 // large enough anyway.
192 if (size == 0)
193 return;
194
195 // If size>0 but c does not have the appropriate
196 // size we throw an exception.
197 //
198 // We could perhaps also allow c.size()>size.
199 // But then looping the loop below get's complicated:
200 // We're not allowed to loop until c.size(). But
201 // we also cannot use size for termination,
202 // because this fails if c is a static vector.
203 if (c.size() != size)
204 DUNE_THROW(RangeError, "Can't resize non-resizable entry v[" << prefix << "] of size " << c.size() << " to size(" << prefix << ")=" << size);
205
206 // Recursively resize all entries of c now.
207 using namespace Dune::Hybrid;
208 prefix.push_back(0);
209 forEach(integralRange(Hybrid::size(c)), [&](auto&& i) {
210 prefix.back() = i;
211 // Here we'd simply like to call resize(c[i], sizeProvider, prefix);
212 // but even gcc-7 does not except this bus reports
213 // "error: ‘this’ was not captured for this lambda function"
214 // although there's no 'this' because we're in a static method.
215 // Bypassing this by and additional method that does perfect
216 // forwarding allows to workaround this.
217 ISTLVectorBackend<V>::forwardToResize(c[i], sizeProvider, prefix);
218 });
219 }
220
221 template<class C, class SizeProvider,
222 std::enable_if_t<not hasResizeMethod<C>::value, int> = 0,
223 std::enable_if_t<isScalar<C>::value, int> = 0>
224 static void resize(C&& c, const SizeProvider& sizeProvider, typename SizeProvider::SizePrefix prefix)
225 {
226 auto size = sizeProvider.size(prefix);
227 if (size != 0)
228 DUNE_THROW(RangeError, "Can't resize scalar vector entry v[" << prefix << "] to size(" << prefix << ")=" << size);
229 }
230
231 template<class C, class T,
232 std::enable_if_t<std::is_assignable<C&,T>::value, int> = 0>
233 void recursiveAssign(C& c, const T& t)
234 {
235 c = t;
236 }
237
238 template<class C, class T,
239 std::enable_if_t<not std::is_assignable<C&,T>::value, int> = 0>
240 void recursiveAssign(C& c, const T& t)
241 {
242 for(auto&& ci: c)
243 recursiveAssign(ci, t);
244 }
245
246public:
247
248 using Vector = V;
249
250 ISTLVectorBackend(Vector& vector) :
251 vector_(&vector)
252 {}
253
254 template<class SizeProvider>
255 void resize(const SizeProvider& sizeProvider)
256 {
257 auto prefix = typename SizeProvider::SizePrefix();
258 prefix.resize(0);
259 resize(*vector_, sizeProvider, prefix);
260 }
261
262 template<class MultiIndex>
263 decltype(auto) operator[](const MultiIndex& index) const
264 {
265 return resolveDynamicMultiIndex(*vector_, index);
266 }
267
268 template<class MultiIndex>
269 decltype(auto) operator[](const MultiIndex& index)
270 {
271 return resolveDynamicMultiIndex(*vector_, index);
272 }
273
282 template<typename T>
283 void operator= (const T& other)
284 {
285 recursiveAssign(vector(), other);
286 }
287
288 template<typename T>
289 void operator= (const ISTLVectorBackend<T>& other)
290 {
291 vector() = other.vector();
292 }
293
294 const Vector& vector() const
295 {
296 return *vector_;
297 }
298
299 Vector& vector()
300 {
301 return *vector_;
302 }
303
304private:
305
306 Vector* vector_;
307};
308
309} // end namespace Impl
310
311
312
344template<class Vector>
345auto istlVectorBackend(Vector& v)
346{
347 static_assert(hasUniqueFieldType<Vector&>(), "Vector type passed to istlVectorBackend() does not have a unique field type.");
348 return Impl::ISTLVectorBackend<Vector>(v);
349}
350
351
352
382template<class Vector>
383auto istlVectorBackend(const Vector& v)
384{
385 static_assert(hasUniqueFieldType<const Vector&>(), "Vector type passed to istlVectorBackend() does not have a unique field type.");
386 return Impl::ISTLVectorBackend<const Vector>(v);
387}
388
389
390
391} // namespace Functions
392} // namespace Dune
393
394
395#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_ISTLVECTORBACKEND_HH
Default exception class for range errors.
Definition: exceptions.hh:252
Infrastructure for concepts.
decltype(auto) unpackIntegerSequence(F &&f, std::integer_sequence< I, i... > sequence)
Unpack an std::integer_sequence<I,i...> to std::integral_constant<I,i>...
Definition: indices.hh:123
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:51
typename detected_or< nonesuch, Op, Args... >::value_t is_detected
Detects whether Op<Args...> is valid.
Definition: type_traits.hh:326
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
auto istlVectorBackend(const Vector &v)
Return a vector backend wrapping const ISTL like containers.
Definition: istlvectorbackend.hh:383
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:266
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:182
constexpr decltype(auto) resolveDynamicMultiIndex(C &&c, const MultiIndex &multiIndex, const IsFinal &isFinal)
Provide multi-index access by chaining operator[].
Definition: indexaccess.hh:354
Dune namespace.
Definition: alignedallocator.hh:11
constexpr auto uniqueTypeList(TypeList< T... > list)
Remove duplicates from a Dune::TypeList.
Definition: typelist.hh:234
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)