DUNE PDELab (git)

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