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