DUNE PDELab (git)

vectorhelpers.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_PDELAB_BACKEND_ISTL_VECTORHELPERS_HH
4#define DUNE_PDELAB_BACKEND_ISTL_VECTORHELPERS_HH
5
7
9
10#include <dune/pdelab/backend/istl/tags.hh>
11#include <dune/pdelab/backend/istl/descriptors.hh>
12#include <dune/pdelab/finiteelementmap/utility.hh>
13#include <dune/pdelab/ordering/orderingbase.hh>
14
15namespace Dune {
16 namespace PDELab {
17
18 // Recursive accessors for vector entries using tag dispatch
19
20#ifndef DOXYGEN // All of the following functions are mere implementation details
21
22 namespace ISTL {
23
24 template<typename CI, typename Block>
25 typename Block::field_type&
26 access_vector_element(tags::field_vector_1, Block& b, const CI& ci, int i)
27 {
28 // usually we are at the end of the multi-index (-1),
29 // but we might be in a PowerFunctionSpace of size 1,
30 // then we are at the lowest multi-index component (0)
31 assert(i == -1 || i == 0);
32 return b[0];
33 }
34
35 template<typename CI, typename Block>
36 typename Block::field_type&
37 access_vector_element(tags::field_vector_n, Block& b, const CI& ci, int i)
38 {
39 assert(i == 0);
40 return b[ci[0]];
41 }
42
43 template<typename CI, typename Block>
44 typename Block::field_type&
45 access_vector_element(tags::block_vector, Block& b, const CI& ci, int i)
46 {
47 return access_vector_element(container_tag(b[ci[i]]),b[ci[i]],ci,i-1);
48 }
49
50
51 template<typename CI, typename Block>
52 const typename Block::field_type&
53 access_vector_element(tags::field_vector_1, const Block& b, const CI& ci, int i)
54 {
55 // usually we are at the end of the multi-index (-1),
56 // but we might be in a PowerFunctionSpace of size 1,
57 // then we are at the lowest multi-index component (0)
58 assert(i == -1 || i == 0);
59 return b[0];
60 }
61
62 template<typename CI, typename Block>
63 const typename Block::field_type&
64 access_vector_element(tags::field_vector_n, const Block& b, const CI& ci, int i)
65 {
66 assert(i == 0);
67 return b[ci[0]];
68 }
69
70 template<typename CI, typename Block>
71 const typename Block::field_type&
72 access_vector_element(tags::block_vector, const Block& b, const CI& ci, int i)
73 {
74 return access_vector_element(container_tag(b[ci[i]]),b[ci[i]],ci,i-1);
75 }
76
77
78 template<typename Vector>
79 [[deprecated]]
80 void resize_vector(tags::block_vector, Vector& v, std::size_t size, bool copy_values)
81 {
82 v.resize(size);
83 }
84
85 template<typename Vector>
86 [[deprecated]]
87 void resize_vector(tags::field_vector, Vector& v, std::size_t size, bool copy_values)
88 {
89 }
90
91 template<typename DI, typename CI, typename Container>
92 [[deprecated]]
93 void allocate_vector(tags::field_vector, const OrderingBase<DI,CI>& ordering, Container& c)
94 {
95 }
96
97 template<typename DI, typename CI, typename Container>
98 [[deprecated]]
99 void allocate_vector(tags::block_vector, const OrderingBase<DI,CI>& ordering, Container& c)
100 {
101 for (std::size_t i = 0; i < ordering.childOrderingCount(); ++i)
102 {
103 if (ordering.containerBlocked())
104 {
105 resize_vector(container_tag(c[i]),c[i],ordering.childOrdering(i).blockCount(),false);
106 allocate_vector(container_tag(c[i]),ordering.childOrdering(i),c[i]);
107 }
108 else
109 allocate_vector(container_tag(c),ordering.childOrdering(i),c);
110 }
111 }
112
113 template<typename Ordering, typename Container>
114 [[deprecated]]
115 void dispatch_vector_allocation(const Ordering& ordering, Container& c, HierarchicContainerAllocationTag tag)
116 {
117 allocate_vector(container_tag(c),ordering,c);
118 }
119
120 template<typename Ordering, typename Container>
121 [[deprecated]]
122 void dispatch_vector_allocation(const Ordering& ordering, Container& c, FlatContainerAllocationTag tag)
123 {
124 resize_vector(container_tag(c),c,ordering.blockCount(),false);
125 }
126
127
128 // ********************************************************************************
129 // TMPs for deducing ISTL block structure from GFS backends
130 // ********************************************************************************
131
132 // tag dispatch switch on GFS tag for per-node functor - general version
133 template<typename E,typename Node, typename Tag, bool isLeafTag = std::is_base_of<LeafGridFunctionSpaceTag,Tag>::value >
134 struct vector_descriptor_helper
135 {
136 // export backend type, as the actual TMP is in the parent reduction functor
137 typedef Node type;
138 };
139
140 // descriptor for backends of leaf spaces collecting various information about
141 // possible blocking structures
142 template<typename E, typename GFS>
143 struct leaf_vector_descriptor
144 {
145
146 using Backend = typename GFS::Traits::Backend;
147 using FEM = typename GFS::Traits::FiniteElementMap;
148
149 static_assert(Backend::Traits::block_type != Blocking::bcrs,
150 "Dynamically blocked leaf spaces are not supported by this backend.");
151
152 // flag for sibling reduction - always true in the leaf case
153 static const bool support_no_blocking = true;
154
155 // flag indicating whether the associated vector type supports cascading
156 // the static blocking further up the tree (i.e. create larger static blocks
157 // at the parent node level. Due to ISTL limitations, this only works once in
158 // the hierarchy, so we only support cascading if we don't already do static
159 // blocking at the current level.
160 static const bool support_cascaded_blocking =
161 Backend::Traits::block_type == Blocking::none; // FIXME
162
163 // The cumulative block size is used by the algorithm to calculate total block
164 // size over several children for cascaded blocking. We try to extract this size
165 // from the finite element map (that works if the number of DOFs per entity is
166 // an identical constant for all types of entities that have DOFs attached), and
167 // if we fail we fall back to 1, which always works.
168 static const std::size_t detected_cumulative_block_size = Dune::Std::detected_or_t<
169 std::integral_constant<std::size_t,0>,
170 FiniteElementMapBlockSize,
171 FEM
172 >::value;
173
174 // Has the user selected an explicit block size? If yes, give priority to that.
175 static const std::size_t cumulative_block_size = Backend::Traits::block_size > 0
176 ? Backend::Traits::block_size
177 : detected_cumulative_block_size;
178
179 static constexpr bool have_valid_block_size = cumulative_block_size > 0;
180
181 static_assert(
182 Backend::Traits::block_size == 0 or (detected_cumulative_block_size % cumulative_block_size) == 0,
183 "The vector block size you specified is not compatible with the finite element map"
184 );
185
186 static_assert(
187 Backend::Traits::block_type != Blocking::fixed or have_valid_block_size,
188 "You requested static blocking, but we cannot extract a valid block size from the finite element map. Please specify the block size with the second template parameter of the vector backend."
189 );
190
191 // The static block size of the associated vector
192 static const std::size_t block_size =
193 Backend::Traits::block_type == Blocking::fixed ? cumulative_block_size : 1;
194
195 // The element type for the vector.
196 typedef E element_type;
197
198 // The ISTL vector type associated with the current subtree.
200
201 };
202
203 // Tag dispatch for leaf spaces - extract leaf descriptor.
204 template<typename E, typename Node, typename Tag>
205 struct vector_descriptor_helper<E,Node,Tag, /* is LeafTag */ true>
206 {
207 typedef leaf_vector_descriptor<E,Node> type;
208 };
209
210 // the actual functor
211 template<typename E>
212 struct extract_vector_descriptor
213 {
214
215 template<typename Node, typename TreePath>
216 struct doVisit
217 {
218 // visit all nodes
219 static const bool value = true;
220 };
221
222 template<typename Node, typename TreePath>
223 struct visit
224 {
225 // forward to actual implementation via tag dispatch
226 typedef typename vector_descriptor_helper<E,Node,TypeTree::ImplementationTag<Node>>::type type;
227 };
228
229 };
230
231 // Descriptor for combining sibling nodes in the tree
232 template<typename Sibling, typename Child>
233 struct cascading_vector_descriptor
234 {
235
236 // We only support cascaded blocking if all children support it
237 static const bool support_cascaded_blocking =
238 Sibling::support_cascaded_blocking &&
239 Child::support_cascaded_blocking;
240
241 // ISTL requires a single, globally constant blocking structure
242 // for its containers, so we make sure the siblings don't disagree
243 // on it.
244 static const bool support_no_blocking =
245 (Sibling::support_no_blocking &&
246 std::is_same<
247 typename Sibling::vector_type,
248 typename Child::vector_type
249 >::value);
250
251 static constexpr bool have_valid_block_size =
252 Sibling::have_valid_block_size and Child::have_valid_block_size;
253
254 // block size
255 static const std::size_t block_size =
256 support_no_blocking ? Sibling::block_size : 1;
257
258 // The element type for the vector.
259 typedef typename Sibling::element_type element_type;
260
261 // Accumulate total block size of all siblings
262 static const std::size_t cumulative_block_size =
263 Sibling::cumulative_block_size + Child::cumulative_block_size;
264
265 // The ISTL vector type associated with the current subtree.
267
268 };
269
270
271 // Switch that turns off standard reduction for the first child of a node.
272 // Default case: do the standard reduction.
273 template<typename D1, typename D2>
274 struct initial_reduction_switch
275 {
276 typedef cascading_vector_descriptor<D1,D2> type;
277 };
278
279 // specialization for first child
280 template<typename D2>
281 struct initial_reduction_switch<void,D2>
282 {
283 typedef D2 type;
284 };
285
286 // sibling reduction functor
287 struct combine_vector_descriptor_siblings
288 {
289
290 template<typename D1, typename D2>
291 struct reduce
292 : public initial_reduction_switch<D1,D2>
293 {};
294
295 };
296
297 // Data part of child -> parent reduction descriptor
298 template<typename Child, typename GFS>
299 struct parent_child_vector_descriptor_data
300 {
301
302 using Backend = typename GFS::Traits::Backend;
303
304 static constexpr bool have_valid_block_size = Child::have_valid_block_size;
305
306 // If all our have a common blocking structure, we can just
307 // concatenate them without doing any blocking
308 static const bool support_no_blocking =
309 Child::support_no_blocking;
310
311 // We support cascaded blocking if neither we nor any of our
312 // children are blocked yet.
313 static const bool support_cascaded_blocking =
314 Child::support_cascaded_blocking &&
315 Backend::Traits::block_type == Blocking::none;
316
317 // It is not allowed to specify a block size on an interior node
318 static_assert(
319 Backend::Traits::block_size == 0,
320 "You cannot specify a block size on interior nodes of the function space tree."
321 );
322
323 // Throw an assertion if the user requests static blocking at this level,
324 // but we cannot support it.
325 static_assert((Backend::Traits::block_type != Blocking::fixed) ||
326 Child::support_cascaded_blocking,
327 "invalid blocking structure.");
328
329 static_assert(
330 Backend::Traits::block_type != Blocking::fixed or have_valid_block_size,
331 "You requested static blocking, but at least one leaf space has a finite element that does not support automatic block size extraction. Please specify the block size with the second template parameter of that space's vector backend."
332 );
333
334 // If we block statically, we create bigger blocks, otherwise the
335 // block size doesn't change.
336 static const std::size_t block_size =
337 Backend::Traits::block_type == Blocking::fixed
338 ? Child::cumulative_block_size
339 : Child::block_size;
340
341 // Just forward this...
342 static const std::size_t cumulative_block_size =
343 Child::cumulative_block_size;
344
345 // The element type for the vector.
346 typedef typename Child::element_type element_type;
347
348 // The ISTL vector type associated with our subtrees.
349 typedef typename Child::vector_type child_vector_type;
350
351 };
352
353 // dispatch switch on blocking type - prototype
354 template<typename Data, Blocking>
355 struct parent_child_vector_descriptor;
356
357 // dispatch switch on blocking type - no blocking case
358 template<typename Data>
359 struct parent_child_vector_descriptor<
360 Data,
361 Blocking::none
362 >
363 : public Data
364 {
365 static_assert(Data::support_no_blocking,
366 "Cannot combine incompatible child block structures without static blocking. "
367 "Did you want to apply static blocking at this level?");
368
369 // Just forward the child vector type
370 typedef typename Data::child_vector_type vector_type;
371 };
372
373 // dispatch switch on blocking type - dynamic blocking case
374 template<typename Data>
375 struct parent_child_vector_descriptor<
376 Data,
377 Blocking::bcrs
378 >
379 : public Data
380 {
381 static_assert(Data::support_no_blocking,
382 "Incompatible child block structures detected, cannot perform dynamic blocking. "
383 "Did you want to apply static blocking at this level?");
384
385 // Wrap the child vector type in another BlockVector
387 };
388
389 // dispatch switch on blocking type - static blocking case
390 template<typename Data>
391 struct parent_child_vector_descriptor<
392 Data,
393 Blocking::fixed
394 >
395 : public Data
396 {
397 // build new block vector with large field block size
398 typedef Dune::BlockVector<
399 FieldVector<
400 typename Data::element_type,
401 Data::block_size
402 >
403 > vector_type;
404 };
405
406 // Child - parent reduction functor
407 struct combine_vector_descriptor_parent
408 {
409
410 template<typename Child, typename GFS>
411 struct reduce
412 {
413
414 struct type
415 : public parent_child_vector_descriptor<parent_child_vector_descriptor_data<
416 Child,
417 GFS>,
418 GFS::Traits::Backend::Traits::block_type
419 >
420 {};
421 };
422
423 };
424
425 // policy describing the GFS tree -> ISTL vector reduction
426 template<typename E>
427 struct vector_creation_policy
428 : public TypeTree::TypeAccumulationPolicy<extract_vector_descriptor<E>,
429 combine_vector_descriptor_siblings,
430 void,
431 combine_vector_descriptor_parent,
432 TypeTree::bottom_up_reduction>
433 {};
434
435 } // namespace ISTL
436
437
438#endif // DOXYGEN
439
440 } // namespace PDELab
441} // namespace Dune
442
443#endif // DUNE_PDELAB_BACKEND_ISTL_VECTORHELPERS_HH
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A vector of blocks with memory management.
Definition: bvector.hh:392
Traits for type conversions and type information.
typename detected_or< Default, Op, Args... >::type detected_or_t
Returns Op<Args...> if that is valid; otherwise returns the fallback type Default.
Definition: type_traits.hh:189
constexpr GeometryType none(unsigned int dim)
Returns a GeometryType representing a singular of dimension dim.
Definition: type.hh:471
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)