dune-functions  2.8.0
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 #include <dune/common/indices.hh>
9 
16 
17 
18 
19 namespace Dune {
20 namespace Functions {
21 
22 
23 // *****************************************************************************
24 // This is the reusable part of the power bases. It contains
25 //
26 // PowerPreBasis
27 //
28 // The pre-basis allows to create the others and is the owner of possible shared
29 // state. These components do _not_ depend on the global basis and local view
30 // and can be used without a global basis.
31 // *****************************************************************************
32 
44 template<class MI, class IMS, class SPB, std::size_t C>
46 {
47  static const std::size_t children = C;
48 
49 public:
50 
52  using SubPreBasis = SPB;
53 
55  using GridView = typename SPB::GridView;
56 
58  using size_type = std::size_t;
59 
61  using IndexMergingStrategy = IMS;
62 
63  using SubNode = typename SubPreBasis::Node;
64 
67 
69  using IndexSet = Impl::DefaultNodeIndexSet<PowerPreBasis>;
70 
72  using MultiIndex = MI;
73 
75  using SizePrefix = Dune::ReservedVector<size_type, MultiIndex::max_size()>;
76 
77 private:
78 
79  using SubMultiIndex = MI;
80 
81 public:
82 
88  template<class... SFArgs,
89  disableCopyMove<PowerPreBasis, SFArgs...> = 0,
90  enableIfConstructible<SubPreBasis, SFArgs...> = 0>
91  PowerPreBasis(SFArgs&&... sfArgs) :
92  subPreBasis_(std::forward<SFArgs>(sfArgs)...)
93  {
94  static_assert(models<Concept::PreBasis<GridView>, SubPreBasis>(), "Subprebasis passed to PowerPreBasis does not model the PreBasis concept.");
95  }
96 
99  {
100  subPreBasis_.initializeIndices();
101  }
102 
104  const GridView& gridView() const
105  {
106  return subPreBasis_.gridView();
107  }
108 
110  void update(const GridView& gv)
111  {
112  subPreBasis_.update(gv);
113  }
114 
118  Node makeNode() const
119  {
120  auto node = Node{};
121  for (std::size_t i=0; i<children; ++i)
122  node.setChild(i, subPreBasis_.makeNode());
123  return node;
124  }
125 
133  [[deprecated("Warning: The IndexSet typedef and the makeIndexSet method are deprecated. "\
134  "As a replacement use the indices() method of the PreBasis directly.")]]
136  {
137  return IndexSet{*this};
138  }
139 
141  size_type size() const
142  {
143  return size({});
144  }
145 
147  size_type size(const SizePrefix& prefix) const
148  {
149  return size(prefix, IndexMergingStrategy{});
150  }
151 
152 private:
153 
155  {
156  // The root index size is the root index size of a single subnode
157  // multiplied by the number of subnodes because we enumerate all
158  // child indices in a row.
159  if (prefix.size() == 0)
160  return children*subPreBasis_.size({});
161 
162  // The first prefix entry refers to one of the (root index size)
163  // subindex trees. Hence we have to first compute the corresponding
164  // prefix entry for a single subnode subnode. The we can append
165  // the other prefix entries unmodified, because the index tree
166  // looks the same after the first level.
167  typename SubPreBasis::SizePrefix subPrefix;
168  subPrefix.push_back(prefix[0] / children);
169  for(std::size_t i=1; i<prefix.size(); ++i)
170  subPrefix.push_back(prefix[i]);
171  return subPreBasis_.size(subPrefix);
172  }
173 
174  size_type size(const SizePrefix& prefix, BasisFactory::FlatLexicographic) const
175  {
176  // The size at the index tree root is the size of at the index tree
177  // root of a single subnode multiplied by the number of subnodes
178  // because we enumerate all child indices in a row.
179  if (prefix.size() == 0)
180  return children*subPreBasis_.size({});
181 
182  // The first prefix entry refers to one of the (root index size)
183  // subindex trees. Hence we have to first compute the corresponding
184  // prefix entry for a single subnode subnode. The we can append
185  // the other prefix entries unmodified, because the index tree
186  // looks the same after the first level.
187  typename SubPreBasis::SizePrefix subPrefix;
188  subPrefix.push_back(prefix[0] % children);
189  for(std::size_t i=1; i<prefix.size(); ++i)
190  subPrefix.push_back(prefix[i]);
191  return subPreBasis_.size(subPrefix);
192  }
193 
194  size_type size(const SizePrefix& prefix, BasisFactory::BlockedLexicographic) const
195  {
196  if (prefix.size() == 0)
197  return children;
198  typename SubPreBasis::SizePrefix subPrefix;
199  for(std::size_t i=1; i<prefix.size(); ++i)
200  subPrefix.push_back(prefix[i]);
201  return subPreBasis_.size(subPrefix);
202  }
203 
204  size_type size(const SizePrefix& prefix, BasisFactory::BlockedInterleaved) const
205  {
206  if (prefix.size() == 0)
207  return subPreBasis_.size();
208 
209  typename SubPreBasis::SizePrefix subPrefix;
210  for(std::size_t i=0; i<prefix.size()-1; ++i)
211  subPrefix.push_back(prefix[i]);
212 
213  size_type r = subPreBasis_.size(subPrefix);
214  if (r==0)
215  return 0;
216  subPrefix.push_back(prefix.back());
217  r = subPreBasis_.size(subPrefix);
218  if (r==0)
219  return children;
220  return r;
221  }
222 
223 public:
224 
227  {
228  return subPreBasis_.dimension() * children;
229  }
230 
233  {
234  return subPreBasis_.maxNodeSize() * children;
235  }
236 
238  const SubPreBasis& subPreBasis() const
239  {
240  return subPreBasis_;
241  }
242 
245  {
246  return subPreBasis_;
247  }
248 
250  template<typename It>
251  It indices(const Node& node, It it) const
252  {
253  return indices(node, it, IndexMergingStrategy{});
254  }
255 
256 private:
257 
258  template<typename It>
259  It indices(const Node& node, It multiIndices, BasisFactory::FlatInterleaved) const
260  {
261  using namespace Dune::Indices;
262  size_type subTreeSize = node.child(_0).size();
263  // Fill indices for first child at the beginning.
264  auto next = Impl::preBasisIndices(subPreBasis(), node.child(_0), multiIndices);
265  // Multiply first component of all indices for first child by
266  // number of children to strech the index range for interleaving.
267  for (std::size_t i = 0; i<subTreeSize; ++i)
268  multiIndices[i][0] *= children;
269  for (std::size_t child = 1; child<children; ++child)
270  {
271  for (std::size_t i = 0; i<subTreeSize; ++i)
272  {
273  // Copy indices from first child for all other children
274  // and shift them by child index to interleave indices.
275  // multiIndices[child*subTreeSize+i] = multiIndices[i];
276  // multiIndices[child*subTreeSize+i][0] = multiIndices[i][0]+child;
277  (*next) = multiIndices[i];
278  (*next)[0] = multiIndices[i][0]+child;
279  ++next;
280  }
281  }
282  return next;
283  }
284 
285  template<typename It>
286  It indices(const Node& node, It multiIndices, BasisFactory::FlatLexicographic) const
287  {
288  using namespace Dune::Indices;
289  size_type subTreeSize = node.child(_0).size();
290  size_type firstIndexEntrySize = subPreBasis().size({});
291  // Fill indices for first child at the beginning.
292  auto next = Impl::preBasisIndices(subPreBasis(), node.child(_0), multiIndices);
293  for (std::size_t child = 1; child<children; ++child)
294  {
295  for (std::size_t i = 0; i<subTreeSize; ++i)
296  {
297  // Copy indices from first child for all other children
298  // and shift them by suitable offset to get lexicographic indices.
299  // multiIndices[child*subTreeSize+i] = multiIndices[i];
300  // multiIndices[child*subTreeSize+i][0] += child*firstIndexEntrySize;
301  (*next) = multiIndices[i];
302  (*next)[0] += child*firstIndexEntrySize;
303  ++next;
304  }
305  }
306  return next;
307  }
308 
309  static void multiIndexPushFront(MultiIndex& M, size_type M0)
310  {
311  M.resize(M.size()+1);
312  for(std::size_t i=M.size()-1; i>0; --i)
313  M[i] = M[i-1];
314  M[0] = M0;
315  }
316 
317  template<typename It>
318  It indices(const Node& node, It multiIndices, BasisFactory::BlockedLexicographic) const
319  {
320  using namespace Dune::Indices;
321  size_type subTreeSize = node.child(_0).size();
322  // Fill indices for first child at the beginning.
323  auto next = Impl::preBasisIndices(subPreBasis(), node.child(_0), multiIndices);
324  // Insert 0 before first component of all indices for first child.
325  for (std::size_t i = 0; i<subTreeSize; ++i)
326  multiIndexPushFront(multiIndices[i], 0);
327  for (std::size_t child = 1; child<children; ++child)
328  {
329  for (std::size_t i = 0; i<subTreeSize; ++i)
330  {
331  // Copy indices from first child for all other children and overwrite
332  // zero in first component as inserted above by child index.
333  // multiIndices[child*subTreeSize+i] = multiIndices[i];
334  // multiIndices[child*subTreeSize+i][0] = child;
335  (*next) = multiIndices[i];
336  (*next)[0] = child;
337  ++next;
338  }
339  }
340  return next;
341  }
342 
343  template<typename It>
344  It indices(const Node& node, It multiIndices, BasisFactory::BlockedInterleaved) const
345  {
346  using namespace Dune::Indices;
347  size_type subTreeSize = node.child(_0).size();
348  // Fill indices for first child at the beginning.
349  auto next = Impl::preBasisIndices(subPreBasis(), node.child(_0), multiIndices);
350  // Append 0 after last component of all indices for first child.
351  for (std::size_t i = 0; i<subTreeSize; ++i)
352  multiIndices[i].push_back(0);
353  for (std::size_t child = 1; child<children; ++child)
354  {
355  for (std::size_t i = 0; i<subTreeSize; ++i)
356  {
357  // Copy indices from first child for all other children and overwrite
358  // zero in last component as appended above by child index.
359  (*next) = multiIndices[i];
360  (*next).back() = child;
361  ++next;
362  }
363  }
364  return next;
365  }
366 
367  SubPreBasis subPreBasis_;
368 };
369 
370 
371 
372 namespace BasisFactory {
373 
374 namespace Imp {
375 
376 template<std::size_t k, class IndexMergingStrategy, class ChildPreBasisFactory>
377 class PowerPreBasisFactory
378 {
379  static const bool isBlocked = std::is_same<IndexMergingStrategy,BlockedLexicographic>::value or std::is_same<IndexMergingStrategy,BlockedInterleaved>::value;
380 
381  static const std::size_t maxChildIndexSize = ChildPreBasisFactory::requiredMultiIndexSize;
382 
383 public:
384 
385  static const std::size_t requiredMultiIndexSize = isBlocked ? (maxChildIndexSize+1) : maxChildIndexSize;
386 
387  PowerPreBasisFactory(const ChildPreBasisFactory& childPreBasisFactory) :
388  childPreBasisFactory_(childPreBasisFactory)
389  {}
390 
391  PowerPreBasisFactory(ChildPreBasisFactory&& childPreBasisFactory) :
392  childPreBasisFactory_(std::move(childPreBasisFactory))
393  {}
394 
395  template<class MultiIndex, class GridView>
396  auto makePreBasis(const GridView& gridView) const
397  {
398  auto childPreBasis = childPreBasisFactory_.template makePreBasis<MultiIndex>(gridView);
399  using ChildPreBasis = decltype(childPreBasis);
400 
401  return PowerPreBasis<MultiIndex, IndexMergingStrategy, ChildPreBasis, k>(std::move(childPreBasis));
402  }
403 
404 private:
405  ChildPreBasisFactory childPreBasisFactory_;
406 };
407 
408 } // end namespace BasisFactory::Imp
409 
410 
411 
424 template<std::size_t k, class ChildPreBasisFactory, class IndexMergingStrategy>
425 auto power(ChildPreBasisFactory&& childPreBasisFactory, const IndexMergingStrategy& ims)
426 {
427  return Imp::PowerPreBasisFactory<k, IndexMergingStrategy, ChildPreBasisFactory>(std::forward<ChildPreBasisFactory>(childPreBasisFactory));
428 }
429 
440 template<std::size_t k, class ChildPreBasisFactory>
441 auto power(ChildPreBasisFactory&& childPreBasisFactory)
442 {
443  return Imp::PowerPreBasisFactory<k, BlockedInterleaved, ChildPreBasisFactory>(std::forward<ChildPreBasisFactory>(childPreBasisFactory));
444 }
445 
446 } // end namespace BasisFactory
447 
448 // Backward compatibility
449 namespace BasisBuilder {
450 
451  using namespace BasisFactory;
452 
453 }
454 
455 
456 } // end namespace Functions
457 } // end namespace Dune
458 
459 
460 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_POWERBASIS_HH
auto power(ChildPreBasisFactory &&childPreBasisFactory)
Create a factory builder that can build a PowerPreBasis.
Definition: powerbasis.hh:441
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:182
Definition: nodes.hh:191
A pre-basis for power bases.
Definition: powerbasis.hh:46
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: powerbasis.hh:232
const SubPreBasis & subPreBasis() const
Const access to the stored prebasis of the factor in the power space.
Definition: powerbasis.hh:238
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: powerbasis.hh:110
void initializeIndices()
Initialize the global indices.
Definition: powerbasis.hh:98
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: powerbasis.hh:104
std::size_t size_type
Type used for indices and size information.
Definition: powerbasis.hh:58
typename SubPreBasis::Node SubNode
Definition: powerbasis.hh:63
size_type size() const
Same as size(prefix) with empty prefix.
Definition: powerbasis.hh:141
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: powerbasis.hh:72
Dune::ReservedVector< size_type, MultiIndex::max_size()> SizePrefix
Type used for prefixes handed to the size() method.
Definition: powerbasis.hh:75
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: powerbasis.hh:226
IMS IndexMergingStrategy
Strategy used to merge the global indices of the child factories.
Definition: powerbasis.hh:61
size_type size(const SizePrefix &prefix) const
Return number of possible values for next position in multi index.
Definition: powerbasis.hh:147
typename SPB::GridView GridView
The grid view that the FE basis is defined on.
Definition: powerbasis.hh:55
SubPreBasis & subPreBasis()
Mutable access to the stored prebasis of the factor in the power space.
Definition: powerbasis.hh:244
Node makeNode() const
Create tree node.
Definition: powerbasis.hh:118
SPB SubPreBasis
The child pre-basis.
Definition: powerbasis.hh:52
PowerPreBasis(SFArgs &&... sfArgs)
Constructor for given child pre-basis objects.
Definition: powerbasis.hh:91
IndexSet makeIndexSet() const
Create tree node index set.
Definition: powerbasis.hh:135
It indices(const Node &node, It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis.
Definition: powerbasis.hh:251
PowerBasisNode< SubNode, children > Node
Template mapping root tree path to type of created tree node.
Definition: powerbasis.hh:66
Impl::DefaultNodeIndexSet< PowerPreBasis > IndexSet
Type of created tree node index set.
Definition: powerbasis.hh:69