dune-functions  2.7.1
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 PB, class IMS>
36 class PowerNodeIndexSet;
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>
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  using SubNode = typename SubPreBasis::Node;
74 
75  using SubIndexSet = typename SubPreBasis::IndexSet;
76 
79 
82 
84  using MultiIndex = MI;
85 
87  using SizePrefix = Dune::ReservedVector<size_type, MultiIndex::max_size()>;
88 
89 private:
90 
91  using SubMultiIndex = MI;
92 
93 public:
94 
100  template<class... SFArgs,
101  disableCopyMove<PowerPreBasis, SFArgs...> = 0,
102  enableIfConstructible<SubPreBasis, SFArgs...> = 0>
103  PowerPreBasis(SFArgs&&... sfArgs) :
104  subPreBasis_(std::forward<SFArgs>(sfArgs)...)
105  {
106  static_assert(models<Concept::PreBasis<GridView>, SubPreBasis>(), "Subprebasis passed to PowerPreBasis does not model the PreBasis concept.");
107  }
108 
111  {
112  subPreBasis_.initializeIndices();
113  }
114 
116  const GridView& gridView() const
117  {
118  return subPreBasis_.gridView();
119  }
120 
122  void update(const GridView& gv)
123  {
124  subPreBasis_.update(gv);
125  }
126 
130  Node makeNode() const
131  {
132  auto node = Node{};
133  for (std::size_t i=0; i<children; ++i)
134  node.setChild(i, subPreBasis_.makeNode());
135  return node;
136  }
137 
145  {
146  return IndexSet{*this};
147  }
148 
150  size_type size() const
151  {
152  return size({});
153  }
154 
156  size_type size(const SizePrefix& prefix) const
157  {
158  return size(prefix, IndexMergingStrategy{});
159  }
160 
161 private:
162 
164  {
165  // The root index size is the root index size of a single subnode
166  // multiplied by the number of subnodes because we enumerate all
167  // child indices in a row.
168  if (prefix.size() == 0)
169  return children*subPreBasis_.size({});
170 
171  // The first prefix entry refers to one of the (root index size)
172  // subindex trees. Hence we have to first compute the corresponding
173  // prefix entry for a single subnode subnode. The we can append
174  // the other prefix entries unmodified, because the index tree
175  // looks the same after the first level.
176  typename SubPreBasis::SizePrefix subPrefix;
177  subPrefix.push_back(prefix[0] / children);
178  for(std::size_t i=1; i<prefix.size(); ++i)
179  subPrefix.push_back(prefix[i]);
180  return subPreBasis_.size(subPrefix);
181  }
182 
183  size_type size(const SizePrefix& prefix, BasisFactory::FlatLexicographic) const
184  {
185  // The size at the index tree root is the size of at the index tree
186  // root of a single subnode multiplied by the number of subnodes
187  // because we enumerate all child indices in a row.
188  if (prefix.size() == 0)
189  return children*subPreBasis_.size({});
190 
191  // The first prefix entry refers to one of the (root index size)
192  // subindex trees. Hence we have to first compute the corresponding
193  // prefix entry for a single subnode subnode. The we can append
194  // the other prefix entries unmodified, because the index tree
195  // looks the same after the first level.
196  typename SubPreBasis::SizePrefix subPrefix;
197  subPrefix.push_back(prefix[0] % children);
198  for(std::size_t i=1; i<prefix.size(); ++i)
199  subPrefix.push_back(prefix[i]);
200  return subPreBasis_.size(subPrefix);
201  }
202 
203  size_type size(const SizePrefix& prefix, BasisFactory::BlockedLexicographic) const
204  {
205  if (prefix.size() == 0)
206  return children;
207  typename SubPreBasis::SizePrefix subPrefix;
208  for(std::size_t i=1; i<prefix.size(); ++i)
209  subPrefix.push_back(prefix[i]);
210  return subPreBasis_.size(subPrefix);
211  }
212 
213  size_type size(const SizePrefix& prefix, BasisFactory::BlockedInterleaved) const
214  {
215  if (prefix.size() == 0)
216  return subPreBasis_.size();
217 
218  typename SubPreBasis::SizePrefix subPrefix;
219  for(std::size_t i=0; i<prefix.size()-1; ++i)
220  subPrefix.push_back(prefix[i]);
221 
222  size_type r = subPreBasis_.size(subPrefix);
223  if (r==0)
224  return 0;
225  subPrefix.push_back(prefix.back());
226  r = subPreBasis_.size(subPrefix);
227  if (r==0)
228  return children;
229  return r;
230  }
231 
232 public:
233 
236  {
237  return subPreBasis_.dimension() * children;
238  }
239 
242  {
243  return subPreBasis_.maxNodeSize() * children;
244  }
245 
247  const SubPreBasis& subPreBasis() const
248  {
249  return subPreBasis_;
250  }
251 
254  {
255  return subPreBasis_;
256  }
257 
258 private:
259  SubPreBasis subPreBasis_;
260 };
261 
262 
263 
264 template<class PB, class IMS>
266 {
267 public:
268 
269  using size_type = std::size_t;
270  using PreBasis = PB;
271  using MultiIndex = typename PreBasis::MultiIndex;
272  using Node = typename PreBasis::Node;
273 
274 protected:
275 
276  using IndexMergingStrategy = IMS;
277  using SubIndexSet = typename PreBasis::SubPreBasis::IndexSet;
278  static const std::size_t children = PreBasis::children;
279 
280 public:
281 
282  PowerNodeIndexSet(const PreBasis & preBasis) :
283  preBasis_(&preBasis),
284  subNodeIndexSet_(preBasis_->subPreBasis().makeIndexSet())
285  {}
286 
287  void bind(const Node& node)
288  {
289  using namespace TypeTree::Indices;
290  node_ = &node;
291  subNodeIndexSet_.bind(node.child(_0));
292  }
293 
294  void unbind()
295  {
296  node_ = nullptr;
297  subNodeIndexSet_.unbind();
298  }
299 
300  size_type size() const
301  {
302  return node_->size();
303  }
304 
306  template<typename It>
307  It indices(It it) const
308  {
309  return indices(it, IndexMergingStrategy{});
310  }
311 
312 private:
313 
314  template<typename It>
315  It indices(It multiIndices, BasisFactory::FlatInterleaved) const
316  {
317  using namespace Dune::TypeTree::Indices;
318  size_type subTreeSize = node_->child(_0).size();
319  // Fill indices for first child at the beginning.
320  auto next = subNodeIndexSet_.indices(multiIndices);
321  // Multiply first component of all indices for first child by
322  // number of children to strech the index range for interleaving.
323  for (std::size_t i = 0; i<subTreeSize; ++i)
324  multiIndices[i][0] *= children;
325  for (std::size_t child = 1; child<children; ++child)
326  {
327  for (std::size_t i = 0; i<subTreeSize; ++i)
328  {
329  // Copy indices from first child for all other children
330  // and shift them by child index to interleave indices.
331  // multiIndices[child*subTreeSize+i] = multiIndices[i];
332  // multiIndices[child*subTreeSize+i][0] = multiIndices[i][0]+child;
333  (*next) = multiIndices[i];
334  (*next)[0] = multiIndices[i][0]+child;
335  ++next;
336  }
337  }
338  return next;
339  }
340 
341  template<typename It>
342  It indices(It multiIndices, BasisFactory::FlatLexicographic) const
343  {
344  using namespace Dune::TypeTree::Indices;
345  size_type subTreeSize = node_->child(_0).size();
346  size_type firstIndexEntrySize = preBasis_->subPreBasis().size({});
347  // Fill indices for first child at the beginning.
348  auto next = subNodeIndexSet_.indices(multiIndices);
349  for (std::size_t child = 1; child<children; ++child)
350  {
351  for (std::size_t i = 0; i<subTreeSize; ++i)
352  {
353  // Copy indices from first child for all other children
354  // and shift them by suitable offset to get lexicographic indices.
355  // multiIndices[child*subTreeSize+i] = multiIndices[i];
356  // multiIndices[child*subTreeSize+i][0] += child*firstIndexEntrySize;
357  (*next) = multiIndices[i];
358  (*next)[0] += child*firstIndexEntrySize;
359  ++next;
360  }
361  }
362  return next;
363  }
364 
365  static void multiIndexPushFront(MultiIndex& M, size_type M0)
366  {
367  M.resize(M.size()+1);
368  for(std::size_t i=M.size()-1; i>0; --i)
369  M[i] = M[i-1];
370  M[0] = M0;
371  }
372 
373  template<typename It>
374  It indices(It multiIndices, BasisFactory::BlockedLexicographic) const
375  {
376  using namespace Dune::TypeTree::Indices;
377  size_type subTreeSize = node_->child(_0).size();
378  // Fill indices for first child at the beginning.
379  auto next = subNodeIndexSet_.indices(multiIndices);
380  // Insert 0 before first component of all indices for first child.
381  for (std::size_t i = 0; i<subTreeSize; ++i)
382  multiIndexPushFront(multiIndices[i], 0);
383  for (std::size_t child = 1; child<children; ++child)
384  {
385  for (std::size_t i = 0; i<subTreeSize; ++i)
386  {
387  // Copy indices from first child for all other children and overwrite
388  // zero in first component as inserted above by child index.
389  // multiIndices[child*subTreeSize+i] = multiIndices[i];
390  // multiIndices[child*subTreeSize+i][0] = child;
391  (*next) = multiIndices[i];
392  (*next)[0] = child;
393  ++next;
394  }
395  }
396  return next;
397  }
398 
399  template<typename It>
400  It indices(It multiIndices, BasisFactory::BlockedInterleaved) const
401  {
402  using namespace Dune::TypeTree::Indices;
403  size_type subTreeSize = node_->child(_0).size();
404  // Fill indices for first child at the beginning.
405  auto next = subNodeIndexSet_.indices(multiIndices);
406  // Append 0 after last component of all indices for first child.
407  for (std::size_t i = 0; i<subTreeSize; ++i)
408  multiIndices[i].push_back(0);
409  for (std::size_t child = 1; child<children; ++child)
410  {
411  for (std::size_t i = 0; i<subTreeSize; ++i)
412  {
413  // Copy indices from first child for all other children and overwrite
414  // zero in last component as appended above by child index.
415  (*next) = multiIndices[i];
416  (*next).back() = child;
417  ++next;
418  }
419  }
420  return next;
421  }
422 
423  const PreBasis* preBasis_;
424  SubIndexSet subNodeIndexSet_;
425  const Node* node_;
426 };
427 
428 
429 
430 namespace BasisFactory {
431 
432 namespace Imp {
433 
434 template<std::size_t k, class IndexMergingStrategy, class ChildPreBasisFactory>
435 class PowerPreBasisFactory
436 {
437  static const bool isBlocked = std::is_same<IndexMergingStrategy,BlockedLexicographic>::value or std::is_same<IndexMergingStrategy,BlockedInterleaved>::value;
438 
439  static const std::size_t maxChildIndexSize = ChildPreBasisFactory::requiredMultiIndexSize;
440 
441 public:
442 
443  static const std::size_t requiredMultiIndexSize = isBlocked ? (maxChildIndexSize+1) : maxChildIndexSize;
444 
445  PowerPreBasisFactory(const ChildPreBasisFactory& childPreBasisFactory) :
446  childPreBasisFactory_(childPreBasisFactory)
447  {}
448 
449  PowerPreBasisFactory(ChildPreBasisFactory&& childPreBasisFactory) :
450  childPreBasisFactory_(std::move(childPreBasisFactory))
451  {}
452 
453  template<class MultiIndex, class GridView>
454  auto makePreBasis(const GridView& gridView) const
455  {
456  auto childPreBasis = childPreBasisFactory_.template makePreBasis<MultiIndex>(gridView);
457  using ChildPreBasis = decltype(childPreBasis);
458 
459  return PowerPreBasis<MultiIndex, IndexMergingStrategy, ChildPreBasis, k>(std::move(childPreBasis));
460  }
461 
462 private:
463  ChildPreBasisFactory childPreBasisFactory_;
464 };
465 
466 } // end namespace BasisFactory::Imp
467 
468 
469 
482 template<std::size_t k, class ChildPreBasisFactory, class IndexMergingStrategy>
483 auto power(ChildPreBasisFactory&& childPreBasisFactory, const IndexMergingStrategy& ims)
484 {
485  return Imp::PowerPreBasisFactory<k, IndexMergingStrategy, ChildPreBasisFactory>(std::forward<ChildPreBasisFactory>(childPreBasisFactory));
486 }
487 
498 template<std::size_t k, class ChildPreBasisFactory>
499 auto power(ChildPreBasisFactory&& childPreBasisFactory)
500 {
501  return Imp::PowerPreBasisFactory<k, BlockedInterleaved, ChildPreBasisFactory>(std::forward<ChildPreBasisFactory>(childPreBasisFactory));
502 }
503 
504 } // end namespace BasisFactory
505 
506 // Backward compatibility
507 namespace BasisBuilder {
508 
509  using namespace BasisFactory;
510 
511 }
512 
513 
514 } // end namespace Functions
515 } // end namespace Dune
516 
517 
518 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_POWERBASIS_HH
auto power(ChildPreBasisFactory &&childPreBasisFactory)
Create a factory builder that can build a PowerPreBasis.
Definition: powerbasis.hh:499
typename std::enable_if< std::is_constructible< T, Args... >::value, int >::type enableIfConstructible
Helper to constrain forwarding constructors.
Definition: type_traits.hh:26
Definition: polynomial.hh:10
Base class for index merging strategies to simplify detection.
Definition: basistags.hh:44
Interleaved merging of direct children without blocking.
Definition: basistags.hh:114
Definition: functionspacebases/concepts.hh:144
Definition: nodes.hh:189
Definition: powerbasis.hh:266
std::size_t size_type
Definition: powerbasis.hh:269
IMS IndexMergingStrategy
Definition: powerbasis.hh:276
PowerNodeIndexSet(const PreBasis &preBasis)
Definition: powerbasis.hh:282
PB PreBasis
Definition: powerbasis.hh:270
size_type size() const
Definition: powerbasis.hh:300
void unbind()
Definition: powerbasis.hh:294
typename PreBasis::MultiIndex MultiIndex
Definition: powerbasis.hh:271
typename PreBasis::Node Node
Definition: powerbasis.hh:272
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:307
void bind(const Node &node)
Definition: powerbasis.hh:287
typename PreBasis::SubPreBasis::IndexSet SubIndexSet
Definition: powerbasis.hh:277
A pre-basis for power bases.
Definition: powerbasis.hh:53
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: powerbasis.hh:241
const SubPreBasis & subPreBasis() const
Const access to the stored prebasis of the factor in the power space.
Definition: powerbasis.hh:247
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: powerbasis.hh:122
void initializeIndices()
Initialize the global indices.
Definition: powerbasis.hh:110
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: powerbasis.hh:116
std::size_t size_type
Type used for indices and size information.
Definition: powerbasis.hh:68
typename SubPreBasis::Node SubNode
Definition: powerbasis.hh:73
size_type size() const
Same as size(prefix) with empty prefix.
Definition: powerbasis.hh:150
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: powerbasis.hh:84
Dune::ReservedVector< size_type, MultiIndex::max_size()> SizePrefix
Type used for prefixes handed to the size() method.
Definition: powerbasis.hh:87
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: powerbasis.hh:235
IMS IndexMergingStrategy
Strategy used to merge the global indices of the child factories.
Definition: powerbasis.hh:71
typename SubPreBasis::IndexSet SubIndexSet
Definition: powerbasis.hh:75
size_type size(const SizePrefix &prefix) const
Return number of possible values for next position in multi index.
Definition: powerbasis.hh:156
typename SPB::GridView GridView
The grid view that the FE basis is defined on.
Definition: powerbasis.hh:65
SubPreBasis & subPreBasis()
Mutable access to the stored prebasis of the factor in the power space.
Definition: powerbasis.hh:253
Node makeNode() const
Create tree node.
Definition: powerbasis.hh:130
SPB SubPreBasis
The child pre-basis.
Definition: powerbasis.hh:62
PowerPreBasis(SFArgs &&... sfArgs)
Constructor for given child pre-basis objects.
Definition: powerbasis.hh:103
IndexSet makeIndexSet() const
Create tree node index set.
Definition: powerbasis.hh:144