dune-pdelab  2.5-dev
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  void resize_vector(tags::block_vector, Vector& v, std::size_t size, bool copy_values)
80  {
81  v.resize(size,copy_values);
82  }
83 
84  template<typename Vector>
85  void resize_vector(tags::field_vector, Vector& v, std::size_t size, bool copy_values)
86  {
87  }
88 
89  template<typename DI, typename CI, typename Container>
90  void allocate_vector(tags::field_vector, const OrderingBase<DI,CI>& ordering, Container& c)
91  {
92  }
93 
94  template<typename DI, typename CI, typename Container>
95  void allocate_vector(tags::block_vector, const OrderingBase<DI,CI>& ordering, Container& c)
96  {
97  for (std::size_t i = 0; i < ordering.childOrderingCount(); ++i)
98  {
99  if (ordering.containerBlocked())
100  {
101  resize_vector(container_tag(c[i]),c[i],ordering.childOrdering(i).blockCount(),false);
102  allocate_vector(container_tag(c[i]),ordering.childOrdering(i),c[i]);
103  }
104  else
105  allocate_vector(container_tag(c),ordering.childOrdering(i),c);
106  }
107  }
108 
109  template<typename Ordering, typename Container>
110  void dispatch_vector_allocation(const Ordering& ordering, Container& c, HierarchicContainerAllocationTag tag)
111  {
112  allocate_vector(container_tag(c),ordering,c);
113  }
114 
115  template<typename Ordering, typename Container>
116  void dispatch_vector_allocation(const Ordering& ordering, Container& c, FlatContainerAllocationTag tag)
117  {
118  resize_vector(container_tag(c),c,ordering.blockCount(),false);
119  }
120 
121 
122  // ********************************************************************************
123  // TMPs for deducing ISTL block structure from GFS backends
124  // ********************************************************************************
125 
126  // tag dispatch switch on GFS tag for per-node functor - general version
128  struct vector_descriptor_helper
129  {
130  // export backend type, as the actual TMP is in the parent reduction functor
131  typedef Node type;
132  };
133 
134  // descriptor for backends of leaf spaces collecting various information about
135  // possible blocking structures
136  template<typename E, typename GFS>
137  struct leaf_vector_descriptor
138  {
139 
140  using Backend = typename GFS::Traits::Backend;
141  using FEM = typename GFS::Traits::FiniteElementMap;
142 
143  static_assert(Backend::Traits::block_type != Blocking::bcrs,
144  "Dynamically blocked leaf spaces are not supported by this backend.");
145 
146  // flag for sibling reduction - always true in the leaf case
147  static const bool support_no_blocking = true;
148 
149  // flag indicating whether the associated vector type supports cascading
150  // the static blocking further up the tree (i.e. create larger static blocks
151  // at the parent node level. Due to ISTL limitations, this only works once in
152  // the hierarchy, so we only support cascading if we don't already do static
153  // blocking at the current level.
154  static const bool support_cascaded_blocking =
155  Backend::Traits::block_type == Blocking::none; // FIXME
156 
157  // The cumulative block size is used by the algorithm to calculate total block
158  // size over several children for cascaded blocking. We try to extract this size
159  // from the finite element map (that works if the number of DOFs per entity is
160  // an identical constant for all types of entities that have DOFs attached), and
161  // if we fail we fall back to 1, which always works.
162  static const std::size_t detected_cumulative_block_size = Dune::Std::detected_or_t<
163  std::integral_constant<std::size_t,0>,
165  FEM
166  >::value;
167 
168  // Has the user selected an explicit block size? If yes, give priority to that.
169  static const std::size_t cumulative_block_size = Backend::Traits::block_size > 0
170  ? Backend::Traits::block_size
171  : detected_cumulative_block_size;
172 
173  static constexpr bool have_valid_block_size = cumulative_block_size > 0;
174 
175  static_assert(
176  Backend::Traits::block_size == 0 or (detected_cumulative_block_size % cumulative_block_size) == 0,
177  "The vector block size you specified is not compatible with the finite element map"
178  );
179 
180  static_assert(
181  Backend::Traits::block_type != Blocking::fixed or have_valid_block_size,
182  "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."
183  );
184 
185  // The static block size of the associated vector
186  static const std::size_t block_size =
187  Backend::Traits::block_type == Blocking::fixed ? cumulative_block_size : 1;
188 
189  // The element type for the vector.
190  typedef E element_type;
191 
192  // The ISTL vector type associated with the current subtree.
193  typedef Dune::BlockVector<FieldVector<E,block_size> > vector_type;
194 
195  };
196 
197  // Tag dispatch for leaf spaces - extract leaf descriptor.
198  template<typename E, typename Node, typename Tag>
199  struct vector_descriptor_helper<E,Node,Tag, /* is LeafTag */ true>
200  {
201  typedef leaf_vector_descriptor<E,Node> type;
202  };
203 
204  // the actual functor
205  template<typename E>
206  struct extract_vector_descriptor
207  {
208 
209  template<typename Node, typename TreePath>
210  struct doVisit
211  {
212  // visit all nodes
213  static const bool value = true;
214  };
215 
216  template<typename Node, typename TreePath>
217  struct visit
218  {
219  // forward to actual implementation via tag dispatch
220  typedef typename vector_descriptor_helper<E,Node,TypeTree::ImplementationTag<Node>>::type type;
221  };
222 
223  };
224 
225  // Descriptor for combining sibling nodes in the tree
226  template<typename Sibling, typename Child>
227  struct cascading_vector_descriptor
228  {
229 
230  // We only support cascaded blocking if all children support it
231  static const bool support_cascaded_blocking =
232  Sibling::support_cascaded_blocking &&
233  Child::support_cascaded_blocking;
234 
235  // ISTL requires a single, globally constant blocking structure
236  // for its containers, so we make sure the siblings don't disagree
237  // on it.
238  static const bool support_no_blocking =
239  (Sibling::support_no_blocking &&
240  std::is_same<
241  typename Sibling::vector_type,
242  typename Child::vector_type
243  >::value);
244 
245  static constexpr bool have_valid_block_size =
246  Sibling::have_valid_block_size and Child::have_valid_block_size;
247 
248  // block size
249  static const std::size_t block_size =
250  support_no_blocking ? Sibling::block_size : 1;
251 
252  // The element type for the vector.
253  typedef typename Sibling::element_type element_type;
254 
255  // Accumulate total block size of all siblings
256  static const std::size_t cumulative_block_size =
257  Sibling::cumulative_block_size + Child::cumulative_block_size;
258 
259  // The ISTL vector type associated with the current subtree.
260  typedef Dune::BlockVector<FieldVector<element_type,block_size> > vector_type;
261 
262  };
263 
264 
265  // Switch that turns off standard reduction for the first child of a node.
266  // Default case: do the standard reduction.
267  template<typename D1, typename D2>
268  struct initial_reduction_switch
269  {
270  typedef cascading_vector_descriptor<D1,D2> type;
271  };
272 
273  // specialization for first child
274  template<typename D2>
275  struct initial_reduction_switch<void,D2>
276  {
277  typedef D2 type;
278  };
279 
280  // sibling reduction functor
281  struct combine_vector_descriptor_siblings
282  {
283 
284  template<typename D1, typename D2>
285  struct reduce
286  : public initial_reduction_switch<D1,D2>
287  {};
288 
289  };
290 
291  // Data part of child -> parent reduction descriptor
292  template<typename Child, typename GFS>
293  struct parent_child_vector_descriptor_data
294  {
295 
296  using Backend = typename GFS::Traits::Backend;
297 
298  static constexpr bool have_valid_block_size = Child::have_valid_block_size;
299 
300  // If all our have a common blocking structure, we can just
301  // concatenate them without doing any blocking
302  static const bool support_no_blocking =
303  Child::support_no_blocking;
304 
305  // We support cascaded blocking if neither we nor any of our
306  // children are blocked yet.
307  static const bool support_cascaded_blocking =
308  Child::support_cascaded_blocking &&
309  Backend::Traits::block_type == Blocking::none;
310 
311  // It is not allowed to specify a block size on an interior node
312  static_assert(
313  Backend::Traits::block_size == 0,
314  "You cannot specify a block size on interior nodes of the function space tree."
315  );
316 
317  // Throw an assertion if the user requests static blocking at this level,
318  // but we cannot support it.
319  static_assert((Backend::Traits::block_type != Blocking::fixed) ||
320  Child::support_cascaded_blocking,
321  "invalid blocking structure.");
322 
323  static_assert(
324  Backend::Traits::block_type != Blocking::fixed or have_valid_block_size,
325  "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."
326  );
327 
328  // If we block statically, we create bigger blocks, otherwise the
329  // block size doesn't change.
330  static const std::size_t block_size =
331  Backend::Traits::block_type == Blocking::fixed
332  ? Child::cumulative_block_size
333  : Child::block_size;
334 
335  // Just forward this...
336  static const std::size_t cumulative_block_size =
337  Child::cumulative_block_size;
338 
339  // The element type for the vector.
340  typedef typename Child::element_type element_type;
341 
342  // The ISTL vector type associated with our subtrees.
343  typedef typename Child::vector_type child_vector_type;
344 
345  };
346 
347  // dispatch switch on blocking type - prototype
348  template<typename Data, Blocking>
349  struct parent_child_vector_descriptor;
350 
351  // dispatch switch on blocking type - no blocking case
352  template<typename Data>
353  struct parent_child_vector_descriptor<
354  Data,
355  Blocking::none
356  >
357  : public Data
358  {
359  static_assert(Data::support_no_blocking,
360  "Cannot combine incompatible child block structures without static blocking. "
361  "Did you want to apply static blocking at this level?");
362 
363  // Just forward the child vector type
364  typedef typename Data::child_vector_type vector_type;
365  };
366 
367  // dispatch switch on blocking type - dynamic blocking case
368  template<typename Data>
369  struct parent_child_vector_descriptor<
370  Data,
371  Blocking::bcrs
372  >
373  : public Data
374  {
375  static_assert(Data::support_no_blocking,
376  "Incompatible child block structures detected, cannot perform dynamic blocking. "
377  "Did you want to apply static blocking at this level?");
378 
379  // Wrap the child vector type in another BlockVector
380  typedef Dune::BlockVector<typename Data::child_vector_type> vector_type;
381  };
382 
383  // dispatch switch on blocking type - static blocking case
384  template<typename Data>
385  struct parent_child_vector_descriptor<
386  Data,
387  Blocking::fixed
388  >
389  : public Data
390  {
391  // build new block vector with large field block size
392  typedef Dune::BlockVector<
393  FieldVector<
394  typename Data::element_type,
395  Data::block_size
396  >
397  > vector_type;
398  };
399 
400  // Child - parent reduction functor
401  struct combine_vector_descriptor_parent
402  {
403 
404  template<typename Child, typename GFS>
405  struct reduce
406  {
407 
408  struct type
409  : public parent_child_vector_descriptor<parent_child_vector_descriptor_data<
410  Child,
411  GFS>,
412  GFS::Traits::Backend::Traits::block_type
413  >
414  {};
415  };
416 
417  };
418 
419  // policy describing the GFS tree -> ISTL vector reduction
420  template<typename E>
421  struct vector_creation_policy
422  : public TypeTree::TypeAccumulationPolicy<extract_vector_descriptor<E>,
423  combine_vector_descriptor_siblings,
424  void,
425  combine_vector_descriptor_parent,
426  TypeTree::bottom_up_reduction>
427  {};
428 
429  } // namespace ISTL
430 
431 
432 #endif // DOXYGEN
433 
434  } // namespace PDELab
435 } // namespace Dune
436 
437 #endif // DUNE_PDELAB_BACKEND_ISTL_VECTORHELPERS_HH
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:246
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:93
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
Create fixed size blocks that each group together a fixed number of DOFs from each child space...
Creates one macro block for each child space, each block is a BlockVector / BCRS matrix.
Blocking
The type of blocking employed at this node in the function space tree.
Definition: istl/descriptors.hh:26
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
No blocking at this level.