dune-pdelab  2.7-git
vectorhelpers.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_PDELAB_BACKEND_ISTL_VECTORHELPERS_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_VECTORHELPERS_HH
5 
6 #include <dune/common/typetraits.hh>
7 
8 #include <dune/istl/bvector.hh>
9 
14 
15 namespace 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
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>,
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.
199  typedef Dune::BlockVector<FieldVector<E,block_size> > vector_type;
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.
266  typedef Dune::BlockVector<FieldVector<element_type,block_size> > vector_type;
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
386  typedef Dune::BlockVector<typename Data::child_vector_type> vector_type;
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
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
std::integral_constant< std::size_t, finiteElementMapBlockSize< FEM >()> FiniteElementMapBlockSize
An alias template that encapsulates the result of finiteElementMapBlockSize<FEM>() in an integral con...
Definition: finiteelementmap/utility.hh:88
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:234
Blocking
The type of blocking employed at this node in the function space tree.
Definition: istl/descriptors.hh:27
@ none
No blocking at this level.
@ bcrs
Creates one macro block for each child space, each block is a BlockVector / BCRS matrix.
@ fixed
Create fixed size blocks that each group together a fixed number of DOFs from each child space.
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139