dune-functions  2.9.0
istlvectorbackend.hh
Go to the documentation of this file.
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>
13 #include <dune/common/concept.hh>
14 
17 
18 
19 namespace Dune {
20 namespace Functions {
21 
22 namespace Impl {
23 
24 template<class V,
25  std::enable_if_t<not Dune::models<Imp::Concept::HasStaticIndexAccess, V>() , int> = 0>
26 auto fieldTypes(V&& /*v*/)
27 {
28  return TypeList<V>{};
29 }
30 
31 template<class V,
32  std::enable_if_t<Dune::models<Imp::Concept::HasStaticIndexAccess, V>(), int> = 0>
33 auto 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 
62 template<class V>
63 constexpr auto fieldTypes()
64 {
65  return decltype(Impl::fieldTypes(std::declval<V>())){};
66 }
67 
73 template<class V>
74 constexpr bool hasUniqueFieldType()
75 {
76  return std::tuple_size<std::decay_t<decltype(fieldTypes<V>())>>::value==1;
77 }
78 
79 
80 
81 namespace 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  */
100 template<class V>
101 class 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>
118  using hasDynamicIndexAccess = Dune::Std::is_detected<dynamicIndexAccess_t, std::remove_reference_t<C>>;
119 
120  template<class C>
121  using hasStaticIndexAccess = Dune::Std::is_detected<staticIndexAccess_t, std::remove_reference_t<C>>;
122 
123  template<class C>
124  using hasResizeMethod = Dune::Std::is_detected<resizeMethod_t, std::remove_reference_t<C>>;
125 
126  template<class C>
127  using isDynamicVector = Dune::Std::is_detected<dynamicIndexAccess_t, std::remove_reference_t<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 gets 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&&, 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  Dune::Hybrid::forEach(c, [&](auto&& ci) {
243  recursiveAssign(ci, t);
244  });
245  }
246 
247 public:
248 
249  using Vector = V;
250 
251  ISTLVectorBackend(Vector& vector) :
252  vector_(&vector)
253  {}
254 
255  template<class SizeProvider>
256  void resize(const SizeProvider& sizeProvider)
257  {
258  auto prefix = typename SizeProvider::SizePrefix();
259  prefix.resize(0);
260  resize(*vector_, sizeProvider, prefix);
261  }
262 
263  template<class MultiIndex>
264  decltype(auto) operator[](const MultiIndex& index) const
265  {
266  return resolveDynamicMultiIndex(*vector_, index);
267  }
268 
269  template<class MultiIndex>
270  decltype(auto) operator[](const MultiIndex& index)
271  {
272  return resolveDynamicMultiIndex(*vector_, index);
273  }
274 
283  template<typename T>
284  void operator= (const T& other)
285  {
286  recursiveAssign(vector(), other);
287  }
288 
289  template<typename T>
290  void operator= (const ISTLVectorBackend<T>& other)
291  {
292  vector() = other.vector();
293  }
294 
295  const Vector& vector() const
296  {
297  return *vector_;
298  }
299 
300  Vector& vector()
301  {
302  return *vector_;
303  }
304 
305 private:
306 
307  Vector* vector_;
308 };
309 
310 } // end namespace Impl
311 
312 
313 
345 template<class Vector>
346 auto istlVectorBackend(Vector& v)
347 {
348  static_assert(hasUniqueFieldType<Vector&>(), "Vector type passed to istlVectorBackend() does not have a unique field type.");
349  return Impl::ISTLVectorBackend<Vector>(v);
350 }
351 
352 
353 
383 template<class Vector>
384 auto istlVectorBackend(const Vector& v)
385 {
386  static_assert(hasUniqueFieldType<const Vector&>(), "Vector type passed to istlVectorBackend() does not have a unique field type.");
387  return Impl::ISTLVectorBackend<const Vector>(v);
388 }
389 
390 
391 
392 } // namespace Functions
393 } // namespace Dune
394 
395 
396 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_ISTLVECTORBACKEND_HH
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition: istlvectorbackend.hh:346
constexpr decltype(auto) resolveDynamicMultiIndex(C &&c, const MultiIndex &multiIndex, const IsFinal &isFinal)
Provide multi-index access by chaining operator[].
Definition: indexaccess.hh:354
Definition: polynomial.hh:10
constexpr auto fieldTypes()
Generate list of field types in container.
Definition: istlvectorbackend.hh:63
constexpr bool hasUniqueFieldType()
Check if container has a unique field type.
Definition: istlvectorbackend.hh:74