dune-functions  2.6-dev
powerbasis.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_FUNCTIONS_FUNCTIONSPACEBASES_POWERBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_POWERBASIS_HH
5 
6 #include <dune/common/reservedvector.hh>
7 #include <dune/common/typeutilities.hh>
8 
9 #include <dune/typetree/powernode.hh>
10 #include <dune/typetree/utility.hh>
11 
17 
18 
19 
20 namespace Dune {
21 namespace Functions {
22 
23 
24 // *****************************************************************************
25 // This is the reusable part of the power bases. It contains
26 //
27 // PowerPreBasis
28 // PowerNodeIndexSet
29 //
30 // The pre-basis allows to create the others and is the owner of possible shared
31 // state. These components do _not_ depend on the global basis or index
32 // set and can be used without a global basis.
33 // *****************************************************************************
34 
35 template<class MI, class TP, class IMS, class SPB, std::size_t C>
37 
38 
39 
51 template<class MI, class IMS, class SPB, std::size_t C>
53 {
54  static const std::size_t children = C;
55 
56  template<class, class, class, class, std::size_t>
57  friend class PowerNodeIndexSet;
58 
59 public:
60 
62  using SubPreBasis = SPB;
63 
65  using GridView = typename SPB::GridView;
66 
68  using size_type = std::size_t;
69 
71  using IndexMergingStrategy = IMS;
72 
73  template<class TP>
74  using SubNode = typename SubPreBasis::template Node<decltype(TypeTree::push_back(TP(), 0))>;
75 
76  template<class TP>
77  using SubIndexSet = typename SubPreBasis::template IndexSet<decltype(TypeTree::push_back(TP(), 0))>;
78 
80  template<class TP>
82 
84  template<class TP>
86 
88  using MultiIndex = MI;
89 
91  using SizePrefix = Dune::ReservedVector<size_type, SubPreBasis::SizePrefix::max_size()+1>;
92 
93 private:
94 
95  using SubMultiIndex = MI;
96 
97 public:
98 
104  template<class... SFArgs,
105  disableCopyMove<PowerPreBasis, SFArgs...> = 0,
106  enableIfConstructible<SubPreBasis, SFArgs...> = 0>
107  PowerPreBasis(SFArgs&&... sfArgs) :
108  subPreBasis_(std::forward<SFArgs>(sfArgs)...)
109  {
110  static_assert(models<Concept::PreBasis<GridView>, SubPreBasis>(), "Subprebasis passed to PowerPreBasis does not model the PreBasis concept.");
111  }
112 
115  {
116  subPreBasis_.initializeIndices();
117  }
118 
120  const GridView& gridView() const
121  {
122  return subPreBasis_.gridView();
123  }
124 
126  void update(const GridView& gv)
127  {
128  subPreBasis_.update(gv);
129  }
130 
141  template<class TP>
142  Node<TP> node(const TP& tp) const
143  {
144  auto node = Node<TP>(tp);
145  for (std::size_t i=0; i<children; ++i)
146  node.setChild(i, subPreBasis_.node(TypeTree::push_back(tp, i)));
147  return node;
148  }
149 
159  template<class TP>
161  {
162  return IndexSet<TP>{*this};
163  }
164 
166  size_type size() const
167  {
168  return size({});
169  }
170 
172  size_type size(const SizePrefix& prefix) const
173  {
174  return size(prefix, IndexMergingStrategy{});
175  }
176 
177 private:
178 
180  {
181  // The root index size is the root index size of a single subnode
182  // multiplied by the number of subnodes because we enumerate all
183  // child indices in a row.
184  if (prefix.size() == 0)
185  return children*subPreBasis_.size({});
186 
187  // The first prefix entry refers to one of the (root index size)
188  // subindex trees. Hence we have to first compute the corresponding
189  // prefix entry for a single subnode subnode. The we can append
190  // the other prefix entries unmodified, because the index tree
191  // looks the same after the first level.
192  typename SubPreBasis::SizePrefix subPrefix;
193  subPrefix.push_back(prefix[0] / children);
194  for(std::size_t i=1; i<prefix.size(); ++i)
195  subPrefix.push_back(prefix[i]);
196  return subPreBasis_.size(subPrefix);
197  }
198 
199  size_type size(const SizePrefix& prefix, BasisBuilder::FlatLexicographic) const
200  {
201  // The size at the index tree root is the size of at the index tree
202  // root of a single subnode multiplied by the number of subnodes
203  // because we enumerate all child indices in a row.
204  if (prefix.size() == 0)
205  return children*subPreBasis_.size({});
206 
207  // The first prefix entry refers to one of the (root index size)
208  // subindex trees. Hence we have to first compute the corresponding
209  // prefix entry for a single subnode subnode. The we can append
210  // the other prefix entries unmodified, because the index tree
211  // looks the same after the first level.
212  typename SubPreBasis::SizePrefix subPrefix;
213  subPrefix.push_back(prefix[0] % children);
214  for(std::size_t i=1; i<prefix.size(); ++i)
215  subPrefix.push_back(prefix[i]);
216  return subPreBasis_.size(subPrefix);
217  }
218 
219  size_type size(const SizePrefix& prefix, BasisBuilder::BlockedLexicographic) const
220  {
221  if (prefix.size() == 0)
222  return children;
223  typename SubPreBasis::SizePrefix subPrefix;
224  for(std::size_t i=1; i<prefix.size(); ++i)
225  subPrefix.push_back(prefix[i]);
226  return subPreBasis_.size(subPrefix);
227  }
228 
229  size_type size(const SizePrefix& prefix, BasisBuilder::LeafBlockedInterleaved) const
230  {
231  if (prefix.size() == 0)
232  return subPreBasis_.size();
233 
234  typename SubPreBasis::SizePrefix subPrefix;
235  for(std::size_t i=0; i<prefix.size()-1; ++i)
236  subPrefix.push_back(prefix[i]);
237 
238  size_type r = subPreBasis_.size(subPrefix);
239  if (r==0)
240  return 0;
241  subPrefix.push_back(prefix.back());
242  r = subPreBasis_.size(subPrefix);
243  if (r==0)
244  return children;
245  return r;
246  }
247 
248 public:
249 
252  {
253  return subPreBasis_.dimension() * children;
254  }
255 
258  {
259  return subPreBasis_.maxNodeSize() * children;
260  }
261 
262 protected:
264 };
265 
266 
267 
268 template<class MI, class TP, class IMS, class SPB, std::size_t C>
269 class PowerNodeIndexSet
270 {
271  static const std::size_t children = C;
272 
273 public:
274 
275  using SubPreBasis = SPB;
276 
278  using GridView = typename SPB::GridView;
279  using size_type = std::size_t;
280  using IndexMergingStrategy = IMS;
281 
283  using MultiIndex = MI;
284 
286 
287  using Node = typename PreBasis::template Node<TP>;
288 
289  using SubTreePath = typename TypeTree::Child<Node,0>::TreePath;
290 
291  using SubNodeIndexSet = typename PreBasis::SubPreBasis::template IndexSet<SubTreePath>;
292 
293  PowerNodeIndexSet(const PreBasis & preBasis) :
294  preBasis_(&preBasis),
295  subNodeIndexSet_(preBasis_->subPreBasis_.template indexSet<SubTreePath>())
296  {}
297 
298  void bind(const Node& node)
299  {
300  using namespace TypeTree::Indices;
301  node_ = &node;
302  subNodeIndexSet_.bind(node.child(_0));
303  }
304 
305  void unbind()
306  {
307  node_ = nullptr;
308  subNodeIndexSet_.unbind();
309  }
310 
311  size_type size() const
312  {
313  return node_->size();
314  }
315 
317  template<typename It>
318  It indices(It it) const
319  {
320  return indices(it, IndexMergingStrategy{});
321  }
322 
323  template<typename It>
324  It indices(It multiIndices, BasisBuilder::FlatInterleaved) const
325  {
326  using namespace Dune::TypeTree::Indices;
327  size_type subTreeSize = node_->child(_0).size();
328  // Fill indices for first child at the beginning.
329  auto next = subNodeIndexSet_.indices(multiIndices);
330  // Multiply first component of all indices for first child by
331  // number of children to strech the index range for interleaving.
332  for (std::size_t i = 0; i<subTreeSize; ++i)
333  multiIndices[i][0] *= children;
334  for (std::size_t child = 1; child<children; ++child)
335  {
336  for (std::size_t i = 0; i<subTreeSize; ++i)
337  {
338  // Copy indices from first child for all other children
339  // and shift them by child index to interleave indices.
340  // multiIndices[child*subTreeSize+i] = multiIndices[i];
341  // multiIndices[child*subTreeSize+i][0] = multiIndices[i][0]+child;
342  (*next) = multiIndices[i];
343  (*next)[0] = multiIndices[i][0]+child;
344  ++next;
345  }
346  }
347  return next;
348  }
349 
350  template<typename It>
351  It indices(It multiIndices, BasisBuilder::FlatLexicographic) const
352  {
353  using namespace Dune::TypeTree::Indices;
354  size_type subTreeSize = node_->child(_0).size();
355  size_type firstIndexEntrySize = preBasis_->subPreBasis_.size({});
356  // Fill indices for first child at the beginning.
357  auto next = subNodeIndexSet_.indices(multiIndices);
358  for (std::size_t child = 1; child<children; ++child)
359  {
360  for (std::size_t i = 0; i<subTreeSize; ++i)
361  {
362  // Copy indices from first child for all other children
363  // and shift them by suitable offset to get lexicographic indices.
364  // multiIndices[child*subTreeSize+i] = multiIndices[i];
365  // multiIndices[child*subTreeSize+i][0] += child*firstIndexEntrySize;
366  (*next) = multiIndices[i];
367  (*next)[0] += child*firstIndexEntrySize;
368  ++next;
369  }
370  }
371  return next;
372  }
373 
374  static const void multiIndexPushFront(MultiIndex& M, size_type M0)
375  {
376  M.resize(M.size()+1);
377  for(std::size_t i=M.size()-1; i>0; --i)
378  M[i] = M[i-1];
379  M[0] = M0;
380  }
381 
382  template<typename It>
383  It indices(It multiIndices, BasisBuilder::BlockedLexicographic) const
384  {
385  using namespace Dune::TypeTree::Indices;
386  size_type subTreeSize = node_->child(_0).size();
387  // Fill indices for first child at the beginning.
388  auto next = subNodeIndexSet_.indices(multiIndices);
389  // Insert 0 before first component of all indices for first child.
390  for (std::size_t i = 0; i<subTreeSize; ++i)
391  multiIndexPushFront(multiIndices[i], 0);
392  for (std::size_t child = 1; child<children; ++child)
393  {
394  for (std::size_t i = 0; i<subTreeSize; ++i)
395  {
396  // Copy indices from first child for all other children and overwrite
397  // zero in first component as inserted above by child index.
398  // multiIndices[child*subTreeSize+i] = multiIndices[i];
399  // multiIndices[child*subTreeSize+i][0] = child;
400  (*next) = multiIndices[i];
401  (*next)[0] = child;
402  ++next;
403  }
404  }
405  return next;
406  }
407 
408  template<typename It>
410  {
411  using namespace Dune::TypeTree::Indices;
412  size_type subTreeSize = node_->child(_0).size();
413  // Fill indices for first child at the beginning.
414  auto next = subNodeIndexSet_.indices(multiIndices);
415  // Append 0 after last component of all indices for first child.
416  for (std::size_t i = 0; i<subTreeSize; ++i)
417  multiIndices[i].push_back(0);
418  for (std::size_t child = 1; child<children; ++child)
419  {
420  for (std::size_t i = 0; i<subTreeSize; ++i)
421  {
422  // Copy indices from first child for all other children and overwrite
423  // zero in last component as appended above by child index.
424  (*next) = multiIndices[i];
425  (*next).back() = child;
426  ++next;
427  }
428  }
429  return next;
430  }
431 
432 private:
433  const PreBasis* preBasis_;
434  SubNodeIndexSet subNodeIndexSet_;
435  const Node* node_;
436 };
437 
438 
439 
440 namespace BasisBuilder {
441 
442 namespace Imp {
443 
444 template<std::size_t k, class IndexMergingStrategy, class ChildPreBasisFactory>
445 class PowerPreBasisFactory
446 {
447  static const bool isBlocked = std::is_same<IndexMergingStrategy,BlockedLexicographic>::value or std::is_same<IndexMergingStrategy,LeafBlockedInterleaved>::value;
448 
449  static const std::size_t maxChildIndexSize = ChildPreBasisFactory::requiredMultiIndexSize;
450 
451 public:
452 
453  static const std::size_t requiredMultiIndexSize = isBlocked ? (maxChildIndexSize+1) : maxChildIndexSize;
454 
455  PowerPreBasisFactory(const ChildPreBasisFactory& childPreBasisFactory) :
456  childPreBasisFactory_(childPreBasisFactory)
457  {}
458 
459  PowerPreBasisFactory(ChildPreBasisFactory&& childPreBasisFactory) :
460  childPreBasisFactory_(std::move(childPreBasisFactory))
461  {}
462 
463  template<class MultiIndex, class GridView>
464  auto makePreBasis(const GridView& gridView) const
465  {
466  auto childPreBasis = childPreBasisFactory_.template makePreBasis<MultiIndex>(gridView);
467  using ChildPreBasis = decltype(childPreBasis);
468 
469  return PowerPreBasis<MultiIndex, IndexMergingStrategy, ChildPreBasis, k>(std::move(childPreBasis));
470  }
471 
472 private:
473  ChildPreBasisFactory childPreBasisFactory_;
474 };
475 
476 } // end namespace BasisBuilder::Imp
477 
478 
479 
492 template<std::size_t k, class ChildPreBasisFactory, class IndexMergingStrategy>
493 auto power(ChildPreBasisFactory&& childPreBasisFactory, const IndexMergingStrategy& ims)
494 {
495  return Imp::PowerPreBasisFactory<k, IndexMergingStrategy, ChildPreBasisFactory>(std::forward<ChildPreBasisFactory>(childPreBasisFactory));
496 }
497 
508 template<std::size_t k, class ChildPreBasisFactory>
509 auto power(ChildPreBasisFactory&& childPreBasisFactory)
510 {
511  return Imp::PowerPreBasisFactory<k, LeafBlockedInterleaved, ChildPreBasisFactory>(std::forward<ChildPreBasisFactory>(childPreBasisFactory));
512 }
513 
514 } // end namespace BasisBuilder
515 
516 
517 
518 } // end namespace Functions
519 } // end namespace Dune
520 
521 
522 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_POWERBASIS_HH
Dune::Functions::PowerPreBasis::SubNode
typename SubPreBasis::template Node< decltype(TypeTree::push_back(TP(), 0))> SubNode
Definition: powerbasis.hh:74
Dune::Functions::BasisBuilder::IndexMergingStrategy
Base class for index merging strategies to simplify detection.
Definition: basistags.hh:44
Dune::Functions::PowerNodeIndexSet::bind
void bind(const Node &node)
Definition: powerbasis.hh:298
Dune::Functions::BasisBuilder::LeafBlockedInterleaved
Interleaved merging of direct children with blocking (i.e. creating blocks at the leaves containing o...
Definition: basistags.hh:180
Dune
Definition: polynomial.hh:7
Dune::Functions::PowerNodeIndexSet::indices
It indices(It multiIndices, BasisBuilder::FlatInterleaved) const
Definition: powerbasis.hh:324
Dune::Functions::PowerPreBasis::node
Node< TP > node(const TP &tp) const
Create tree node with given root tree path.
Definition: powerbasis.hh:142
Dune::Functions::PowerNodeIndexSet::multiIndexPushFront
static const void multiIndexPushFront(MultiIndex &M, size_type M0)
Definition: powerbasis.hh:374
Dune::Functions::PowerPreBasis::update
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: powerbasis.hh:126
Dune::Functions::PowerNodeIndexSet::unbind
void unbind()
Definition: powerbasis.hh:305
Dune::Functions::PowerNodeIndexSet::SubNodeIndexSet
typename PreBasis::SubPreBasis::template IndexSet< SubTreePath > SubNodeIndexSet
Definition: powerbasis.hh:291
Dune::Functions::PowerNodeIndexSet::IndexMergingStrategy
IMS IndexMergingStrategy
Definition: powerbasis.hh:280
Dune::Functions::PowerNodeIndexSet::indices
It indices(It multiIndices, BasisBuilder::BlockedLexicographic) const
Definition: powerbasis.hh:383
concepts.hh
Dune::Functions::PowerPreBasis::GridView
typename SPB::GridView GridView
The grid view that the FE basis is defined on.
Definition: powerbasis.hh:65
Dune::Functions::PowerPreBasis::size_type
std::size_t size_type
Type used for indices and size information.
Definition: powerbasis.hh:68
Dune::Functions::PowerPreBasis::subPreBasis_
SubPreBasis subPreBasis_
Definition: powerbasis.hh:263
Dune::Functions::BasisBuilder::power
auto power(ChildPreBasisFactory &&childPreBasisFactory)
Create a factory builder that can build a PowerPreBasis.
Definition: powerbasis.hh:509
Dune::Functions::PowerNodeIndexSet::SubPreBasis
SPB SubPreBasis
Definition: powerbasis.hh:275
Dune::Functions::PowerPreBasis::gridView
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: powerbasis.hh:120
Dune::Functions::PowerNodeIndexSet::SubTreePath
typename TypeTree::Child< Node, 0 >::TreePath SubTreePath
Definition: powerbasis.hh:289
Dune::Functions::PowerNodeIndexSet::MultiIndex
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: powerbasis.hh:283
Dune::Functions::PowerPreBasis::SizePrefix
Dune::ReservedVector< size_type, SubPreBasis::SizePrefix::max_size()+1 > SizePrefix
Type used for prefixes handed to the size() method.
Definition: powerbasis.hh:91
Dune::Functions::PowerPreBasis::indexSet
IndexSet< TP > indexSet() const
Create tree node index set with given root tree path.
Definition: powerbasis.hh:160
Dune::Functions::PowerNodeIndexSet::PowerNodeIndexSet
PowerNodeIndexSet(const PreBasis &preBasis)
Definition: powerbasis.hh:293
Dune::Functions::PowerPreBasis::MultiIndex
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: powerbasis.hh:88
nodes.hh
Dune::Functions::PowerNodeIndexSet::indices
It indices(It multiIndices, BasisBuilder::LeafBlockedInterleaved) const
Definition: powerbasis.hh:409
Dune::Functions::PowerNodeIndexSet::size_type
std::size_t size_type
Definition: powerbasis.hh:279
Dune::Functions::PowerPreBasis::initializeIndices
void initializeIndices()
Initialize the global indices.
Definition: powerbasis.hh:114
Dune::Functions::PowerNodeIndexSet::Node
typename PreBasis::template Node< TP > Node
Definition: powerbasis.hh:287
Dune::Functions::PowerPreBasis::dimension
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: powerbasis.hh:251
Dune::Functions::BasisBuilder::FlatLexicographic
Lexicographic merging of direct children without blocking.
Definition: basistags.hh:78
Dune::Functions::PowerNodeIndexSet::indices
It indices(It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis.
Definition: powerbasis.hh:318
Dune::Functions::PowerNodeIndexSet::indices
It indices(It multiIndices, BasisBuilder::FlatLexicographic) const
Definition: powerbasis.hh:351
utility.hh
Dune::Functions::PowerPreBasis::maxNodeSize
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: powerbasis.hh:257
Dune::Functions::PowerPreBasis::PowerPreBasis
PowerPreBasis(SFArgs &&... sfArgs)
Constructor for given child pre-basis objects.
Definition: powerbasis.hh:107
type_traits.hh
Dune::Functions::BasisBuilder::BlockedLexicographic
Lexicographic merging of direct children with blocking (i.e. creating one block per direct child).
Definition: basistags.hh:146
Dune::Functions::PowerPreBasis
A pre-basis for power bases.
Definition: powerbasis.hh:52
Dune::Functions::PowerBasisNode
Definition: nodes.hh:209
Dune::Functions::PowerNodeIndexSet::GridView
typename SPB::GridView GridView
The grid view that the FE space is defined on.
Definition: powerbasis.hh:278
basistags.hh
Dune::Functions::PowerPreBasis::SubPreBasis
SPB SubPreBasis
The child pre-basis.
Definition: powerbasis.hh:62
Dune::Functions::PowerPreBasis::IndexMergingStrategy
IMS IndexMergingStrategy
Strategy used to merge the global indices of the child factories.
Definition: powerbasis.hh:71
Dune::Functions::PowerPreBasis::size
size_type size() const
Same as size(prefix) with empty prefix.
Definition: powerbasis.hh:166
Dune::Functions::Concept::PreBasis
Definition: concepts.hh:153
Dune::Functions::enableIfConstructible
typename std::enable_if< std::is_constructible< T, Args... >::value, int >::type enableIfConstructible
Helper to constrain forwarding constructors.
Definition: type_traits.hh:26
Dune::Functions::BasisBuilder::FlatInterleaved
Interleaved merging of direct children without blocking.
Definition: basistags.hh:112
Dune::Functions::PowerPreBasis::SubIndexSet
typename SubPreBasis::template IndexSet< decltype(TypeTree::push_back(TP(), 0))> SubIndexSet
Definition: powerbasis.hh:77
Dune::Functions::PowerNodeIndexSet
Definition: powerbasis.hh:36
Dune::Functions::PowerNodeIndexSet::size
size_type size() const
Definition: powerbasis.hh:311
Dune::Functions::PowerPreBasis::size
size_type size(const SizePrefix &prefix) const
Return number of possible values for next position in multi index.
Definition: powerbasis.hh:172