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 // this is here for backwards compatibility and deprecation warnings, remove after 2.5.0
7 #include "ensureistlinclude.hh"
8 
9 #include <dune/common/typetraits.hh>
10 
11 #include <dune/istl/bvector.hh>
12 
16 
17 namespace Dune {
18  namespace PDELab {
19 
20  // Recursive accessors for vector entries using tag dispatch
21 
22 #ifndef DOXYGEN // All of the following functions are mere implementation details
23 
24  namespace ISTL {
25 
26  template<typename CI, typename Block>
27  typename Block::field_type&
28  access_vector_element(tags::field_vector_1, Block& b, const CI& ci, int i)
29  {
30  // usually we are at the end of the multi-index (-1),
31  // but we might be in a PowerFunctionSpace of size 1,
32  // then we are at the lowest multi-index component (0)
33  assert(i == -1 || i == 0);
34  return b[0];
35  }
36 
37  template<typename CI, typename Block>
38  typename Block::field_type&
39  access_vector_element(tags::field_vector_n, Block& b, const CI& ci, int i)
40  {
41  assert(i == 0);
42  return b[ci[0]];
43  }
44 
45  template<typename CI, typename Block>
46  typename Block::field_type&
47  access_vector_element(tags::block_vector, Block& b, const CI& ci, int i)
48  {
49  return access_vector_element(container_tag(b[ci[i]]),b[ci[i]],ci,i-1);
50  }
51 
52 
53  template<typename CI, typename Block>
54  const typename Block::field_type&
55  access_vector_element(tags::field_vector_1, const Block& b, const CI& ci, int i)
56  {
57  // usually we are at the end of the multi-index (-1),
58  // but we might be in a PowerFunctionSpace of size 1,
59  // then we are at the lowest multi-index component (0)
60  assert(i == -1 || i == 0);
61  return b[0];
62  }
63 
64  template<typename CI, typename Block>
65  const typename Block::field_type&
66  access_vector_element(tags::field_vector_n, const Block& b, const CI& ci, int i)
67  {
68  assert(i == 0);
69  return b[ci[0]];
70  }
71 
72  template<typename CI, typename Block>
73  const typename Block::field_type&
74  access_vector_element(tags::block_vector, const Block& b, const CI& ci, int i)
75  {
76  return access_vector_element(container_tag(b[ci[i]]),b[ci[i]],ci,i-1);
77  }
78 
79 
80  template<typename Vector>
81  void resize_vector(tags::block_vector, Vector& v, std::size_t size, bool copy_values)
82  {
83  v.resize(size,copy_values);
84  }
85 
86  template<typename Vector>
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  void allocate_vector(tags::field_vector, const OrderingBase<DI,CI>& ordering, Container& c)
93  {
94  }
95 
96  template<typename DI, typename CI, typename Container>
97  void allocate_vector(tags::block_vector, const OrderingBase<DI,CI>& ordering, Container& c)
98  {
99  for (std::size_t i = 0; i < ordering.childOrderingCount(); ++i)
100  {
101  if (ordering.containerBlocked())
102  {
103  resize_vector(container_tag(c[i]),c[i],ordering.childOrdering(i).blockCount(),false);
104  allocate_vector(container_tag(c[i]),ordering.childOrdering(i),c[i]);
105  }
106  else
107  allocate_vector(container_tag(c),ordering.childOrdering(i),c);
108  }
109  }
110 
111  template<typename Ordering, typename Container>
112  void dispatch_vector_allocation(const Ordering& ordering, Container& c, HierarchicContainerAllocationTag tag)
113  {
114  allocate_vector(container_tag(c),ordering,c);
115  }
116 
117  template<typename Ordering, typename Container>
118  void dispatch_vector_allocation(const Ordering& ordering, Container& c, FlatContainerAllocationTag tag)
119  {
120  resize_vector(container_tag(c),c,ordering.blockCount(),false);
121  }
122 
123 
124  // ********************************************************************************
125  // TMPs for deducing ISTL block structure from GFS backends
126  // ********************************************************************************
127 
128  // tag dispatch switch on GFS tag for per-node functor - general version
130  struct vector_descriptor_helper
131  {
132  // export backend type, as the actual TMP is in the parent reduction functor
133  typedef typename Node::Traits::Backend type;
134  };
135 
136  // descriptor for backends of leaf spaces collecting various information about
137  // possible blocking structures
138  template<typename E, typename Backend>
139  struct leaf_vector_descriptor
140  {
141 
142  static_assert(Backend::Traits::block_type != Blocking::bcrs,
143  "Dynamically blocked leaf spaces are not supported by this backend.");
144 
145  // flag for sibling reduction - always true in the leaf case
146  static const bool support_no_blocking = true;
147 
148  // flag indicating whether the associated vector type supports cascading
149  // the static blocking further up the tree (i.e. create larger static blocks
150  // at the parent node level. Due to ISTL limitations, this only works once in
151  // the hierarchy, so we only support cascading if we don't already do static
152  // blocking at the current level.
153  static const bool support_cascaded_blocking =
154  Backend::Traits::block_type == Blocking::none; // FIXME
155 
156  // The static block size of the associated vector
157  static const std::size_t block_size =
158  Backend::Traits::block_type == Blocking::fixed
159  ? Backend::Traits::block_size
160  : 1;
161 
162  // The cumulative block size is used by the algorithm to calculate total block
163  // size over several children for cascaded blocking. Right now, this will always be set to
164  // the block size passed in by the user, but it might also be possible to provide this
165  // information in the FiniteElementMap and provide automatic blocking detection.
166  static const std::size_t cumulative_block_size = Backend::Traits::block_size;
167 
168  // The element type for the vector.
169  typedef E element_type;
170 
171  // The ISTL vector type associated with the current subtree.
172  typedef Dune::BlockVector<FieldVector<E,block_size> > vector_type;
173 
174  };
175 
176  // Tag dispatch for leaf spaces - extract leaf descriptor.
177  template<typename E, typename Node, typename Tag>
178  struct vector_descriptor_helper<E,Node,Tag, /* is LeafTag */ true>
179  {
180  typedef leaf_vector_descriptor<E,typename Node::Traits::Backend> type;
181  };
182 
183  // the actual functor
184  template<typename E>
185  struct extract_vector_descriptor
186  {
187 
188  template<typename Node, typename TreePath>
189  struct doVisit
190  {
191  // visit all nodes
192  static const bool value = true;
193  };
194 
195  template<typename Node, typename TreePath>
196  struct visit
197  {
198  // forward to actual implementation via tag dispatch
199  typedef typename vector_descriptor_helper<E,Node,TypeTree::ImplementationTag<Node>>::type type;
200  };
201 
202  };
203 
204  // Descriptor for combining sibling nodes in the tree
205  template<typename Sibling, typename Child>
206  struct cascading_vector_descriptor
207  {
208 
209  // We only support cascaded blocking if all children support it
210  static const bool support_cascaded_blocking =
211  Sibling::support_cascaded_blocking &&
212  Child::support_cascaded_blocking;
213 
214  // ISTL requires a single, globally constant blocking structure
215  // for its containers, so we make sure the siblings don't disagree
216  // on it.
217  static const bool support_no_blocking =
218  (Sibling::support_no_blocking &&
219  std::is_same<
220  typename Sibling::vector_type,
221  typename Child::vector_type
222  >::value);
223 
224  // block size
225  static const std::size_t block_size =
226  support_no_blocking ? Sibling::block_size : 1;
227 
228  // The element type for the vector.
229  typedef typename Sibling::element_type element_type;
230 
231  // Accumulate total block size of all siblings
232  static const std::size_t cumulative_block_size =
233  Sibling::cumulative_block_size + Child::cumulative_block_size;
234 
235  // The ISTL vector type associated with the current subtree.
236  typedef Dune::BlockVector<FieldVector<element_type,block_size> > vector_type;
237 
238  };
239 
240 
241  // Switch that turns off standard reduction for the first child of a node.
242  // Default case: do the standard reduction.
243  template<typename D1, typename D2>
244  struct initial_reduction_switch
245  {
246  typedef cascading_vector_descriptor<D1,D2> type;
247  };
248 
249  // specialization for first child
250  template<typename D2>
251  struct initial_reduction_switch<void,D2>
252  {
253  typedef D2 type;
254  };
255 
256  // sibling reduction functor
257  struct combine_vector_descriptor_siblings
258  {
259 
260  template<typename D1, typename D2>
261  struct reduce
262  : public initial_reduction_switch<D1,D2>
263  {};
264 
265  };
266 
267  // Data part of child -> parent reduction descriptor
268  template<typename Child, typename Backend>
269  struct parent_child_vector_descriptor_data
270  {
271 
272  // If all our have a common blocking structure, we can just
273  // concatenate them without doing any blocking
274  static const bool support_no_blocking =
275  Child::support_no_blocking;
276 
277  // We support cascaded blocking if neither we nor any of our
278  // children are blocked yet.
279  static const bool support_cascaded_blocking =
280  Child::support_cascaded_blocking &&
281  Backend::Traits::block_type == Blocking::none;
282 
283  // Throw an assertion if the user requests static blocking at this level,
284  // but we cannot support it.
285  static_assert((Backend::Traits::block_type != Blocking::fixed) ||
286  Child::support_cascaded_blocking,
287  "invalid blocking structure.");
288 
289  // If we block statically, we create bigger blocks, otherwise the
290  // block size doesn't change.
291  static const std::size_t block_size =
292  Backend::Traits::block_type == Blocking::fixed
293  ? Child::cumulative_block_size
294  : Child::block_size;
295 
296  // Just forward this...
297  static const std::size_t cumulative_block_size =
298  Child::cumulative_block_size;
299 
300  // The element type for the vector.
301  typedef typename Child::element_type element_type;
302 
303  // The ISTL vector type associated with our subtrees.
304  typedef typename Child::vector_type child_vector_type;
305 
306  };
307 
308  // dispatch switch on blocking type - prototype
309  template<typename Data, Blocking>
310  struct parent_child_vector_descriptor;
311 
312  // dispatch switch on blocking type - no blocking case
313  template<typename Data>
314  struct parent_child_vector_descriptor<
315  Data,
316  Blocking::none
317  >
318  : public Data
319  {
320  static_assert(Data::support_no_blocking,
321  "Cannot combine incompatible child block structures without static blocking. "
322  "Did you want to apply static blocking at this level?");
323 
324  // Just forward the child vector type
325  typedef typename Data::child_vector_type vector_type;
326  };
327 
328  // dispatch switch on blocking type - dynamic blocking case
329  template<typename Data>
330  struct parent_child_vector_descriptor<
331  Data,
332  Blocking::bcrs
333  >
334  : public Data
335  {
336  static_assert(Data::support_no_blocking,
337  "Incompatible child block structures detected, cannot perform dynamic blocking. "
338  "Did you want to apply static blocking at this level?");
339 
340  // Wrap the child vector type in another BlockVector
341  typedef Dune::BlockVector<typename Data::child_vector_type> vector_type;
342  };
343 
344  // dispatch switch on blocking type - static blocking case
345  template<typename Data>
346  struct parent_child_vector_descriptor<
347  Data,
348  Blocking::fixed
349  >
350  : public Data
351  {
352  // build new block vector with large field block size
353  typedef Dune::BlockVector<
354  FieldVector<
355  typename Data::element_type,
356  Data::block_size
357  >
358  > vector_type;
359  };
360 
361  // Child - parent reduction functor
362  struct combine_vector_descriptor_parent
363  {
364 
365  template<typename Child, typename Backend>
366  struct reduce
367  {
368 
369  struct type
370  : public parent_child_vector_descriptor<parent_child_vector_descriptor_data<
371  Child,
372  Backend>,
373  Backend::Traits::block_type
374  >
375  {};
376  };
377 
378  };
379 
380  // policy describing the GFS tree -> ISTL vector reduction
381  template<typename E>
382  struct vector_creation_policy
383  : public TypeTree::TypeAccumulationPolicy<extract_vector_descriptor<E>,
384  combine_vector_descriptor_siblings,
385  void,
386  combine_vector_descriptor_parent,
387  TypeTree::bottom_up_reduction>
388  {};
389 
390  } // namespace ISTL
391 
392 #endif // DOXYGEN
393 
394  } // namespace PDELab
395 } // namespace Dune
396 
397 #endif // DUNE_PDELAB_BACKEND_ISTL_VECTORHELPERS_HH
Create fixed size blocks that each group together a fixed number of DOFs from each child space...
Blocking
The type of blocking employed at this node in the function space tree.
Definition: istl/descriptors.hh:29
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
Creates one macro block for each child space, each block is a BlockVector / BCRS matrix.
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
No blocking at this level.
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:249