dune-functions  2.6-dev
compositebasis.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_COMPOSITEBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_COMPOSITEBASIS_HH
5 
6 #include <tuple>
7 #include <utility>
8 
9 #include <dune/common/std/utility.hh>
10 #include <dune/common/std/apply.hh>
11 #include <dune/common/hybridutilities.hh>
12 #include <dune/common/reservedvector.hh>
13 #include <dune/common/typeutilities.hh>
14 #include <dune/common/hybridutilities.hh>
15 
16 #include <dune/typetree/compositenode.hh>
17 #include <dune/typetree/utility.hh>
18 
25 
26 
27 namespace Dune {
28 namespace Functions {
29 
30 namespace Imp {
31 
32  template<typename... T>
33  struct SizeOf
34  : public std::integral_constant<std::size_t,sizeof...(T)>
35  {};
36 
37  template<typename... T>
38  using index_sequence_for = std::make_index_sequence<SizeOf<T...>{}>;
39 }
40 
41 // *****************************************************************************
42 // This is the reusable part of the composite bases. It contains
43 //
44 // CompositePreBasis
45 // CompositeNodeIndexSet
46 //
47 // The pre-basis allows to create the others and is the owner of possible shared
48 // state. These components do _not_ depend on the global basis or index
49 // set and can be used without a global basis.
50 // *****************************************************************************
51 
52 
53 template<class MI, class TP, class IT, class... SPB>
55 
68 template<class MI, class IMS, class... SPB>
70 {
71 public:
72 
74  using SubPreBases = std::tuple<SPB...>;
75 
77  using GridView = typename std::tuple_element<0, SubPreBases>::type::GridView;
78 
80  using size_type = std::size_t;
81 
83  using IndexMergingStrategy = IMS;
84 
85 protected:
86  static const std::size_t children = sizeof...(SPB);
87 
88  template<class, class, class, class...>
89  friend class CompositeNodeIndexSet;
90 
91  using ChildIndexTuple = IntegerSequenceTuple<Imp::index_sequence_for<SPB...>>;
92 
93  template<class TP>
94  struct FixedTP
95  {
96 
97  template<class I>
98  using IndexToSubTreePath = decltype(TypeTree::push_back(TP(), I()));
99 
101 
102  template<class F, class SubTP>
103  using PreBasisToSubNode = typename F::template Node<SubTP>;
104 
106 
107  template<class F, class SubTP>
108  using PreBasisToSubIndexSet = typename F::template IndexSet<SubTP>;
109 
111 
112  template<class... N>
114 
116  };
117 
118 
119 public:
120 
122  template<std::size_t k>
123  using SubPreBasis = typename std::tuple_element<k, std::tuple<SPB...>>::type;
124 
126  template<class TP>
127  using Node = typename FixedTP<TP>::Node;
128 
130  template<class TP>
131  using IndexSet = CompositeNodeIndexSet<MI, TP, IMS, SPB...>;
132 
134  using MultiIndex = MI;
135 
137  using SizePrefix = Dune::ReservedVector<size_type, MultiIndex::max_size()+1>;
138 
144  template<class... SFArgs,
145  disableCopyMove<CompositePreBasis, SFArgs...> = 0,
146  enableIfConstructible<std::tuple<SPB...>, SFArgs...> = 0>
147  CompositePreBasis(SFArgs&&... sfArgs) :
148  subPreBases_(std::forward<SFArgs>(sfArgs)...)
149  {
150  using namespace Dune::Hybrid;
151  forEach(subPreBases_, [&](const auto& subPreBasis){
152  static_assert(models<Concept::PreBasis<GridView>, std::decay_t<decltype(subPreBasis)>>(), "Subprebases passed to CompositePreBasis does not model the PreBasis concept.");
153  });
154  }
155 
158  {
159  using namespace Dune::Hybrid;
160  forEach(Dune::Std::make_index_sequence<children>(), [&](auto i) {
161  elementAt(subPreBases_, i).initializeIndices();
162  });
163  }
164 
166  const GridView& gridView() const
167  {
168  return std::get<0>(subPreBases_).gridView();
169  }
170 
172  void update(const GridView& gv)
173  {
174  using namespace Dune::Hybrid;
175  forEach(Dune::Std::make_index_sequence<children>(), [&](auto i) {
176  elementAt(subPreBases_, i).update(gv);
177  });
178  }
179 
190  template<class TP>
191  Node<TP> node(const TP& tp) const
192  {
193  auto node = Node<TP>(tp);
194  using namespace Dune::Hybrid;
195  forEach(Dune::Std::make_index_sequence<children>(), [&](auto i) {
196  node.setChild( elementAt(subPreBases_, i).node(TypeTree::push_back(tp, i)), i);
197  });
198  return node;
199  }
200 
210  template<class TP>
212  {
213  return IndexSet<TP>{*this};
214  }
215 
217  size_type size() const
218  {
219  return size({});
220  }
221 
223  size_type size(const SizePrefix& prefix) const
224  {
225  return size(prefix, IndexMergingStrategy{});
226  }
227 
228 private:
229 
230  size_type size(const SizePrefix& prefix, BasisBuilder::BlockedLexicographic) const
231  {
232  if (prefix.size() == 0)
233  return children;
234 
235  return Hybrid::switchCases(std::make_index_sequence<children>(), prefix[0], [&] (auto i) {
236  const auto& subPreBasis = std::get<i.value>(subPreBases_);
237  typename std::decay<decltype(subPreBasis)>::type::SizePrefix subPrefix;
238  for(std::size_t i=1; i<prefix.size(); ++i)
239  subPrefix.push_back(prefix[i]);
240  return subPreBasis.size(subPrefix);
241  }, []() {
242  return size_type(0);
243  });
244  }
245 
246  struct Lambda_size_flat_sizeInSubtree
247  {
248  template<class I, class PB>
249  size_type operator()(const I&, const PB& subPreBases, const SizePrefix& prefix, size_type& shiftedFirst, size_type& r)
250  {
251  using SubPreBasis = typename std::tuple_element<I::value, PB>::type;
252  const SubPreBasis& subPreBasis = std::get<I::value>(subPreBases);
253  if (shiftedFirst < subPreBasis.size())
254  {
255  typename SubPreBasis::SizePrefix subPrefix;
256  subPrefix.push_back(shiftedFirst);
257  for(std::size_t i=1; i<prefix.size(); ++i)
258  subPrefix.push_back(prefix[i]);
259  r = subPreBasis.size(subPrefix);
260  return true;
261  }
262  shiftedFirst -= subPreBasis.size();
263  return false;
264  }
265  };
266 
267  size_type size(const SizePrefix& prefix, BasisBuilder::FlatLexicographic) const
268  {
269  size_type r = 0;
270  using namespace Dune::Hybrid;
271  if (prefix.size() == 0)
272  forEach(Dune::Std::make_index_sequence<children>(), [&](auto i) {
273  r += elementAt(subPreBases_, i).size();
274  });
275  else {
276  size_type shiftedFirst = prefix[0];
277  staticFindInRange<0, sizeof...(SPB)>(Lambda_size_flat_sizeInSubtree(), subPreBases_, prefix, shiftedFirst, r);
278  }
279  return r;
280  }
281 
282 public:
283 
286  {
287  size_type r=0;
288  // Accumulate dimension() for all subprebases
289  using namespace Dune::Hybrid;
290  forEach(Dune::Std::make_index_sequence<children>(), [&](auto i) {
291  r += elementAt(subPreBases_, i).dimension();
292  });
293  return r;
294  }
295 
298  {
299  size_type r=0;
300  // Accumulate maxNodeSize() for all subprebases
301  using namespace Dune::Hybrid;
302  forEach(Dune::Std::make_index_sequence<children>(), [&](auto i) {
303  r += elementAt(subPreBases_, i).maxNodeSize();
304  });
305  return r;
306  }
307 
308 protected:
309  std::tuple<SPB...> subPreBases_;
310 };
311 
312 
313 
314 template<class MI, class TP, class IMS, class... SPB>
316 {
317  static const std::size_t children = sizeof...(SPB);
318 
319 public:
320 
321  template<std::size_t k>
322  using SubPreBasis = typename std::tuple_element<k, std::tuple<SPB...>>::type;
323 
325  using size_type = std::size_t;
326  using IndexMergingStrategy = IMS;
327 
329  using MultiIndex = MI;
330 
331  using PreBasis = CompositePreBasis<MI, IMS, SPB...>;
332 
333  using Node = typename PreBasis::template Node<TP>;
334 
335  using SubTreePaths = typename PreBasis::template FixedTP<TP>::SubTreePaths;
336  using SubIndexSets = typename PreBasis::template FixedTP<TP>::SubIndexSets;
337 
338 
340  {
341  // transform a single (preBasis,subTreePath) pair to subIndexSet
342  template<class SubPreBasis, class SubTP>
343  auto operator()(const SubPreBasis& preBasis, const SubTP& subTP)
344  ->decltype(preBasis.template indexSet<SubTP>())
345  {
346  return preBasis.template indexSet<SubTP>();
347  }
348  };
349 
350  CompositeNodeIndexSet(const PreBasis & preBasis) :
351  preBasis_(&preBasis),
352  subNodeIndexSetTuple_(transformTuple(Lambda_PreBasisToSubIndexSet(), preBasis_->subPreBases_, SubTreePaths()))
353  {}
354 
355  void bind(const Node& node)
356  {
357  node_ = &node;
358  using namespace Dune::Hybrid;
359  forEach(Dune::Std::make_index_sequence<children>(), [&](auto i) {
360  elementAt(subNodeIndexSetTuple_, i).bind(node.child(i));
361  });
362  }
363 
364  void unbind()
365  {
366  node_ = nullptr;
367  using namespace Dune::Hybrid;
368  forEach(Dune::Std::make_index_sequence<children>(), [&](auto i) {
369  elementAt(subNodeIndexSetTuple_, i).unbind();
370  });
371  }
372 
373  size_type size() const
374  {
375  return node_->size();
376  }
377 
379  template<typename It>
380  It indices(It it) const
381  {
382  return indices(it, IndexMergingStrategy{});
383  }
384 
385  template<typename It>
386  It indices(It multiIndices, BasisBuilder::FlatLexicographic) const
387  {
388  using namespace Dune::Hybrid;
389  size_type firstComponentOffset = 0;
390  // Loop over all children
391  forEach(Dune::Std::make_index_sequence<children>(), [&](auto child){
392  const auto& subNodeIndexSet = elementAt(subNodeIndexSetTuple_, child);
393  const auto& subPreBasis = elementAt(preBasis_->subPreBases_, child);
394  size_type subTreeSize = subNodeIndexSet.size();
395  // Fill indices for current child into index buffer starting from current
396  // buffer position and shift first index component of any index for current
397  // child by suitable offset to get lexicographic indices.
398  subNodeIndexSet.indices(multiIndices);
399  for (std::size_t i = 0; i<subTreeSize; ++i)
400  multiIndices[i][0] += firstComponentOffset;
401  // Increment offset by the size for first index component of the current child
402  firstComponentOffset += subPreBasis.size({});
403  // Increment buffer iterator by the number of indices processed for current child
404  multiIndices += subTreeSize;
405  });
406  return multiIndices;
407  }
408 
409  static const void multiIndexPushFront(MultiIndex& M, size_type M0)
410  {
411  M.resize(M.size()+1);
412  for(std::size_t i=M.size()-1; i>0; --i)
413  M[i] = M[i-1];
414  M[0] = M0;
415  }
416 
417  template<typename It>
418  It indices(It multiIndices, BasisBuilder::BlockedLexicographic) const
419  {
420  using namespace Dune::Hybrid;
421  // Loop over all children
422  forEach(Dune::Std::make_index_sequence<children>(), [&](auto child){
423  const auto& subNodeIndexSet = elementAt(subNodeIndexSetTuple_, child);
424  size_type subTreeSize = subNodeIndexSet.size();
425  // Fill indices for current child into index buffer starting from current position
426  subNodeIndexSet.indices(multiIndices);
427  // Insert child index before first component of all indices of current child.
428  for (std::size_t i = 0; i<subTreeSize; ++i)
429  this->multiIndexPushFront(multiIndices[i], child);
430  // Increment buffer iterator by the number of indices processed for current child
431  multiIndices += subTreeSize;
432  });
433  return multiIndices;
434  }
435 
436 
437 private:
438  const PreBasis* preBasis_;
439  SubIndexSets subNodeIndexSetTuple_;
440  const Node* node_;
441 };
442 
443 
444 
445 namespace BasisBuilder {
446 
447 namespace Imp {
448 
449 template<class ST0>
450 constexpr std::size_t maxHelper(ST0&& i0)
451 {
452  return i0;
453 }
454 
455 template<class ST0, class... ST>
456 constexpr std::size_t maxHelper(ST0&& i0, ST&&... i)
457 {
458  return (i0 > maxHelper(i...)) ? i0 : maxHelper(i...);
459 }
460 
461 template<class IndexMergingStrategy, class... ChildPreBasisFactory>
462 class CompositePreBasisFactory
463 {
464  static const bool isBlocked = std::is_same<IndexMergingStrategy,BlockedLexicographic>::value or std::is_same<IndexMergingStrategy,LeafBlockedInterleaved>::value;
465 
466  static const std::size_t maxChildIndexSize = maxHelper(ChildPreBasisFactory::requiredMultiIndexSize...);
467 
468  template<class MultiIndex, class GridView, class... ChildPreBasis>
469  auto makePreBasisFromChildPreBases(const GridView&, ChildPreBasis&&... childPreBasis) const
470  {
471  return CompositePreBasis<MultiIndex, IndexMergingStrategy, std::decay_t<ChildPreBasis>...>(std::forward<ChildPreBasis>(childPreBasis)...);
472  }
473 
474 public:
475 
476  static const std::size_t requiredMultiIndexSize = isBlocked ? (maxChildIndexSize+1) : maxChildIndexSize;
477 
478  CompositePreBasisFactory(const ChildPreBasisFactory&... childPreBasisFactory) :
479  childPreBasisFactories_(childPreBasisFactory...)
480  {}
481 
482  CompositePreBasisFactory(ChildPreBasisFactory&&... childPreBasisFactory) :
483  childPreBasisFactories_(std::move(childPreBasisFactory)...)
484  {}
485 
486  template<class MultiIndex, class GridView>
487  auto makePreBasis(const GridView& gridView) const
488  {
489  // Use Std::apply to unpack the tuple childPreBasisFactories_
490  return Std::apply([&](const auto&... childPreBasisFactory) {
491  return this->makePreBasisFromChildPreBases<MultiIndex>(gridView, childPreBasisFactory.template makePreBasis<MultiIndex>(gridView)...);
492  }, childPreBasisFactories_);
493  }
494 
495 private:
496  std::tuple<ChildPreBasisFactory...> childPreBasisFactories_;
497 };
498 
499 } // end namespace BasisBuilder::Imp
500 
501 
502 
513 template<
514  typename... Args,
515  std::enable_if_t<Concept::isIndexMergingStrategy<typename LastType<Args...>::type>(),int> = 0>
516 auto composite(Args&&... args)
517 {
518  // We have to separate the last entry which is the IndexMergingStrategy
519  // and the preceding ones, which are the ChildPreBasisFactories
520 
521  using ArgTuple = std::tuple<std::decay_t<Args>...>;
522 
523  // Compute number of children and index of the IndexMergingStrategy argument
524  constexpr std::size_t children = Dune::SizeOf<Args...>::value-1;
525 
526  // Use last type as IndexMergingStrategy
527  using IndexMergingStrategy = std::tuple_element_t<children, ArgTuple>;
528 
529  // Index sequence for all but the last entry for partial tuple unpacking
530  auto childIndices = std::make_index_sequence<children>{};
531 
532  // Unpack tuple only for those entries related to children
533  return applyPartial([&](auto&&... childPreBasisFactory){
534  return Imp::CompositePreBasisFactory<IndexMergingStrategy, std::decay_t<decltype(childPreBasisFactory)>...>(std::forward<decltype(childPreBasisFactory)>(childPreBasisFactory)...);
535  },
536  std::forward_as_tuple(std::forward<Args>(args)...),
537  childIndices);
538 }
539 
551 template<
552  typename... Args,
553  std::enable_if_t<not Concept::isIndexMergingStrategy<typename LastType<Args...>::type>(),int> = 0>
554 auto composite(Args&&... args)
555 {
556  return Imp::CompositePreBasisFactory<BasisBuilder::BlockedLexicographic, std::decay_t<Args>...>(std::forward<Args>(args)...);
557 }
558 
559 } // end namespace BasisBuilder
560 
561 
562 
563 } // end namespace Functions
564 } // end namespace Dune
565 
566 
567 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_COMPOSITEBASIS_HH
auto transformTuple(F &&f, const std::tuple< T... > &tuple) -> decltype(Imp::transformTupleHelper(std::forward< F >(f), tuple, std::index_sequence_for< T... >
Transform tuple value using a functor.
Definition: utility.hh:162
TransformTuple< PreBasisToSubNode, SubPreBases, SubTreePaths > SubNodes
Definition: compositebasis.hh:105
void initializeIndices()
Initialize the global indices.
Definition: compositebasis.hh:157
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition: compositebasis.hh:297
ExpandTuple< SubNodesToNode, SubNodes > Node
Definition: compositebasis.hh:115
It indices(It multiIndices, BasisBuilder::FlatLexicographic) const
Definition: compositebasis.hh:386
Definition: compositebasis.hh:94
typename F::template IndexSet< SubTP > PreBasisToSubIndexSet
Definition: compositebasis.hh:108
std::tuple< SPB... > subPreBases_
Definition: compositebasis.hh:309
typename std::tuple_element< k, std::tuple< SPB... > >::type SubPreBasis
Definition: compositebasis.hh:322
static const void multiIndexPushFront(MultiIndex &M, size_type M0)
Definition: compositebasis.hh:409
auto operator()(const SubPreBasis &preBasis, const SubTP &subTP) -> decltype(preBasis.template indexSet< SubTP >())
Definition: compositebasis.hh:343
typename SubPreBasis< 0 >::GridView GridView
Definition: compositebasis.hh:324
void bind(const Node &node)
Definition: compositebasis.hh:355
A pre-basis for composite bases.
Definition: compositebasis.hh:69
auto composite(Args &&... args)
Create a factory builder that can build a CompositePreBasis.
Definition: compositebasis.hh:516
IntegerSequenceTuple< Imp::index_sequence_for< SPB... > > ChildIndexTuple
Definition: compositebasis.hh:91
CompositePreBasis(SFArgs &&... sfArgs)
Constructor for given child pre-basis objects.
Definition: compositebasis.hh:147
IMS IndexMergingStrategy
Strategy used to merge the global indices of the child pre-bases.
Definition: compositebasis.hh:83
Lexicographic merging of direct children without blocking.
Definition: basistags.hh:78
Dune::ReservedVector< size_type, MultiIndex::max_size()+1 > SizePrefix
Type used for prefixes handed to the size() method.
Definition: compositebasis.hh:137
typename F::template Node< SubTP > PreBasisToSubNode
Definition: compositebasis.hh:103
Get last entry of type list.
Definition: utility.hh:218
decltype(TypeTree::push_back(TP(), I())) IndexToSubTreePath
Definition: compositebasis.hh:98
Lexicographic merging of direct children with blocking (i.e. creating one block per direct child)...
Definition: basistags.hh:146
typename FixedTP< TP >::Node Node
Template mapping root tree path to type of created tree node.
Definition: compositebasis.hh:127
Base class for index merging strategies to simplify detection.
Definition: basistags.hh:44
Definition: compositebasis.hh:54
typename Imp::TransformTupleHelper< F, Tuples... >::Type TransformTuple
Transform tuple types argument using type-functor.
Definition: utility.hh:128
std::size_t size_type
Definition: compositebasis.hh:325
void unbind()
Definition: compositebasis.hh:364
CompositeNodeIndexSet(const PreBasis &preBasis)
Definition: compositebasis.hh:350
Definition: nodes.hh:232
TransformTuple< PreBasisToSubIndexSet, SubPreBases, SubTreePaths > SubIndexSets
Definition: compositebasis.hh:110
friend class CompositeNodeIndexSet
Definition: compositebasis.hh:89
Node< TP > node(const TP &tp) const
Create tree node with given root tree path.
Definition: compositebasis.hh:191
Definition: concepts.hh:153
std::size_t size_type
Type used for indices and size information.
Definition: compositebasis.hh:80
decltype(auto) applyPartial(F &&f, ArgTuple &&args, std::integer_sequence< I, i... > indices)
Apply function with arguments from a given tuple.
Definition: utility.hh:267
IndexSet< TP > indexSet() const
Create tree node index set with given root tree path.
Definition: compositebasis.hh:211
typename PreBasis::template FixedTP< TP >::SubIndexSets SubIndexSets
Definition: compositebasis.hh:336
It indices(It multiIndices, BasisBuilder::BlockedLexicographic) const
Definition: compositebasis.hh:418
typename std::tuple_element< k, std::tuple< SPB... > >::type SubPreBasis
Template mapping index of child to its pre-basis type.
Definition: compositebasis.hh:123
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: compositebasis.hh:172
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: compositebasis.hh:166
static constexpr bool isIndexMergingStrategy()
Definition: basistags.hh:23
IMS IndexMergingStrategy
Definition: compositebasis.hh:326
TransformTuple< IndexToSubTreePath, ChildIndexTuple > SubTreePaths
Definition: compositebasis.hh:100
size_type size() const
Definition: compositebasis.hh:373
Definition: polynomial.hh:7
It indices(It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis...
Definition: compositebasis.hh:380
typename Imp::IntegerSequenceTupleHelper< IntegerSequence >::Type IntegerSequenceTuple
Transform integer_sequence<I,k...> to tuple<integral_constant<I,k>...>
Definition: utility.hh:208
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: compositebasis.hh:134
std::tuple< SPB... > SubPreBases
Tuple of child pre-bases.
Definition: compositebasis.hh:74
typename Imp::ExpandTupleHelper< T, ArgTuple >::Type ExpandTuple
Expand tuple arguments as template arguments.
Definition: utility.hh:91
size_type size(const SizePrefix &prefix) const
Return number of possible values for next position in multi index.
Definition: compositebasis.hh:223
size_type size() const
Same as size(prefix) with empty prefix.
Definition: compositebasis.hh:217
typename std::tuple_element< 0, SubPreBases >::type::GridView GridView
The grid view that the FE basis is defined on.
Definition: compositebasis.hh:77
void staticFindInRange(F &&f, Args &&... args)
Static find loop.
Definition: staticforloop.hh:56
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition: compositebasis.hh:285
typename PreBasis::template Node< TP > Node
Definition: compositebasis.hh:333
typename PreBasis::template FixedTP< TP >::SubTreePaths SubTreePaths
Definition: compositebasis.hh:335
STL namespace.
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: compositebasis.hh:329
typename std::enable_if< std::is_constructible< T, Args... >::value, int >::type enableIfConstructible
Helper to constrain forwarding constructors.
Definition: type_traits.hh:26