dune-functions  2.8.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 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 
246 public:
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 
304 private:
305 
306  Vector* vector_;
307 };
308 
309 } // end namespace Impl
310 
311 
312 
344 template<class Vector>
345 auto 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 
382 template<class Vector>
383 auto 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
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition: istlvectorbackend.hh:345
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