DUNE PDELab (2.8)

powerbasis.hh
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
8#include <dune/common/indices.hh>
9
10#include <dune/functions/common/utility.hh>
11#include <dune/functions/common/type_traits.hh>
12#include <dune/functions/functionspacebases/basistags.hh>
13#include <dune/functions/functionspacebases/nodes.hh>
14#include <dune/functions/functionspacebases/concepts.hh>
15#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
16
17
18
19namespace Dune {
20namespace 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
44template<class MI, class IMS, class SPB, std::size_t C>
46{
47 static const std::size_t children = C;
48
49public:
50
52 using SubPreBasis = SPB;
53
55 using GridView = typename SPB::GridView;
56
58 using size_type = std::size_t;
59
62
63 using SubNode = typename SubPreBasis::Node;
64
66 using Node = PowerBasisNode<SubNode, children>;
67
69 using IndexSet = Impl::DefaultNodeIndexSet<PowerPreBasis>;
70
72 using MultiIndex = MI;
73
75 using SizePrefix = Dune::ReservedVector<size_type, MultiIndex::max_size()>;
76
77private:
78
79 using SubMultiIndex = MI;
80
81public:
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
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
142 {
143 return size({});
144 }
145
147 size_type size(const SizePrefix& prefix) const
148 {
149 return size(prefix, IndexMergingStrategy{});
150 }
151
152private:
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
223public:
224
227 {
228 return subPreBasis_.dimension() * children;
229 }
230
233 {
234 return subPreBasis_.maxNodeSize() * children;
235 }
236
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
256private:
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
372namespace BasisFactory {
373
374namespace Imp {
375
376template<std::size_t k, class IndexMergingStrategy, class ChildPreBasisFactory>
377class 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
383public:
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
404private:
405 ChildPreBasisFactory childPreBasisFactory_;
406};
407
408} // end namespace BasisFactory::Imp
409
410
411
424template<std::size_t k, class ChildPreBasisFactory, class IndexMergingStrategy>
425auto power(ChildPreBasisFactory&& childPreBasisFactory, const IndexMergingStrategy& ims)
426{
427 return Imp::PowerPreBasisFactory<k, IndexMergingStrategy, ChildPreBasisFactory>(std::forward<ChildPreBasisFactory>(childPreBasisFactory));
428}
429
440template<std::size_t k, class ChildPreBasisFactory>
441auto 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
449namespace 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
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
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
std::size_t size_type
Type used for indices and size information.
Definition: powerbasis.hh:58
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
SubPreBasis & subPreBasis()
Mutable access to the stored prebasis of the factor in the power space.
Definition: powerbasis.hh:244
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
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: powerbasis.hh:104
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
const SubPreBasis & subPreBasis() const
Const access to the stored prebasis of the factor in the power space.
Definition: powerbasis.hh:238
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
Grid view abstract base class.
Definition: gridview.hh:63
A Vector class with statically reserved memory.
Definition: reservedvector.hh:43
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:51
constexpr HybridTreePath< T..., std::size_t > push_back(const HybridTreePath< T... > &tp, std::size_t i)
Appends a run time index to a HybridTreePath.
Definition: treepath.hh:278
ImplementationDefined child(Node &&node, Indices... indices)
Extracts the child of a node given by a sequence of compile-time and run-time indices.
Definition: childextraction.hh:126
std::enable_if_t< not Impl::disableCopyMoveHelper< This, T... >::value, int > disableCopyMove
Helper to disable constructor as copy and move constructor.
Definition: typeutilities.hh:43
typename std::enable_if< std::is_constructible< T, Args... >::value, int >::type enableIfConstructible
Helper to constrain forwarding constructors.
Definition: type_traits.hh:26
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:49
Dune namespace.
Definition: alignedallocator.hh:11
constexpr Mantissa power(Mantissa m, Exponent p)
Power method for integer exponents.
Definition: math.hh:73
STL namespace.
An stl-compliant random-access container which stores everything on the stack.
Interleaved merging of direct children without blocking.
Definition: basistags.hh:114
Base class for index merging strategies to simplify detection.
Definition: basistags.hh:44
Utilities for type computations, constraining overloads, ...
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)