DUNE PDELab (git)

interleavedordering.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_PDELAB_ORDERING_INTERLEAVEDORDERING_HH
4#define DUNE_PDELAB_ORDERING_INTERLEAVEDORDERING_HH
5
6#include <array>
7#include <string>
8
9#include <dune/typetree/compositenode.hh>
10#include <dune/typetree/powernode.hh>
11
13
14#include <dune/pdelab/gridfunctionspace/tags.hh>
15#include <dune/pdelab/ordering/utility.hh>
16#include <dune/pdelab/ordering/orderingbase.hh>
17
18namespace Dune {
19 namespace PDELab {
20
23
24 namespace interleaved_ordering {
25
27 template<typename DI, typename CI, typename Node>
28 class Base
29 : public OrderingBase<DI,CI>
30 {
31
32 typedef OrderingBase<DI,CI> BaseT;
33
34 public:
35
36 typedef typename OrderingBase<DI,CI>::Traits Traits;
37
39
40 static const bool consume_tree_index = true;
41
43
48 Base(Node& node, bool container_blocked, const OrderingTag& ordering_tag, typename BaseT::GFSData* gfs_data)
49 : BaseT(node,container_blocked,ordering_tag.offsets(),gfs_data,nullptr)
50 {
51 // This check looks a little weird, but there is always one offset more than
52 // there are blocks (the first offsets is 0, and the last one is the "offset
53 // beyond the end" to encode the size of the final child).
54 if (node.degree() + 1 != ordering_tag.offsets().size())
56 "Invalid block structure for InterleavedOrdering: "
57 << node.degree() << " children, but "
58 << (ordering_tag.offsets().size() - 1) << " block sizes.");
59 }
60
61 template<typename ItIn, typename ItOut>
62 void map_lfs_indices(const ItIn begin, const ItIn end, ItOut out) const
63 {
64 typedef typename Traits::SizeType size_type;
65 if (this->_container_blocked)
66 {
67 for (ItIn in = begin; in != end; ++in, ++out)
68 {
69 size_type child_index = in->treeIndex().back();
70 size_type child_block_offset = this->_child_block_merge_offsets[child_index];
71 size_type child_block_size = this->_child_block_merge_offsets[child_index + 1] - child_block_offset;
72 size_type index = out->back();
73 size_type block_index = index / child_block_size;
74 size_type offset = index % child_block_size;
75 size_type block_offset = child_block_offset + offset;
76 out->back() = block_offset;
77 out->push_back(block_index);
78 }
79 }
80 else
81 {
82 for (ItIn in = begin; in != end; ++in, ++out)
83 {
84 size_type child_index = in->treeIndex().back();
85 size_type child_block_offset = this->_child_block_merge_offsets[child_index];
86 size_type child_block_size = this->_child_block_merge_offsets[child_index + 1] - child_block_offset;
87 size_type block_size = this->_child_block_merge_offsets.back();
88 size_type index = out->back();
89 size_type block_index = index / child_block_size;
90 size_type offset = index % child_block_size;
91 size_type block_offset = child_block_offset + offset;
92 out->back() = block_index * block_size + block_offset;
93 }
94 }
95 }
96
97 template<typename CIOutIterator, typename DIOutIterator = DummyDOFIndexIterator>
98 typename Traits::SizeType
99 extract_entity_indices(const typename Traits::DOFIndex::EntityIndex& ei,
100 typename Traits::SizeType child_index,
101 CIOutIterator ci_out, const CIOutIterator ci_end) const
102 {
103 typedef typename Traits::SizeType size_type;
104 if (this->_container_blocked)
105 {
106 for (; ci_out != ci_end; ++ci_out)
107 {
108 size_type child_block_offset = this->_child_block_merge_offsets[child_index];
109 size_type child_block_size = this->_child_block_merge_offsets[child_index + 1] - child_block_offset;
110 size_type index = ci_out->back();
111 size_type block_index = index / child_block_size;
112 size_type offset = index % child_block_size;
113 size_type block_offset = child_block_offset + offset;
114 ci_out->back() = block_offset;
115 ci_out->push_back(block_index);
116 }
117 }
118 else
119 {
120 for (; ci_out != ci_end; ++ci_out)
121 {
122 size_type child_block_offset = this->_child_block_merge_offsets[child_index];
123 size_type child_block_size = this->_child_block_merge_offsets[child_index + 1] - child_block_offset;
124 size_type block_size = this->_child_block_merge_offsets.back();
125 size_type index = ci_out->back();
126 size_type block_index = index / child_block_size;
127 size_type offset =index % child_block_size;
128 ci_out->back() = block_index * block_size + child_block_offset + offset;
129 }
130 }
131
132 // The return value is not used for non-leaf orderings.
133 return 0;
134 }
135
136 using BaseT::size;
137
138 typename Traits::SizeType size(typename Traits::ContainerIndex suffix) const
139 {
140 if (suffix.size() == 0)
141 return this->_block_count;
142
143 // block_size: the interleaving is partitioned into `n` blocks with size `block_size`
144 const auto& block_size = this->_child_block_merge_offsets.back();
145 // block_index: i-th block in this partitioning
146 std::size_t block_index = std::numeric_limits<std::size_t>::max();
147 // block_offset: index within the i-th block
148 std::size_t block_offset = std::numeric_limits<std::size_t>::max();
149
150 if (this->containerBlocked()) {
151 block_index = suffix.back();
152 if (suffix.size() == 1) {
153 assert(block_index < this->_block_count);
154 return block_size;
155 }
156 suffix.pop_back();
157 block_offset = suffix.back();
158 assert(block_offset < block_size);
159 } else {
160 block_index = suffix.back() / block_size;
161 assert(block_index < this->_block_count);
162 block_offset = suffix.back() % block_size;
163 }
164 assert(block_index < this->_child_block_offsets.back());
165
166 const auto& merge_begin = std::begin(this->_child_block_merge_offsets);
167 const auto& merge_end = std::end(this->_child_block_merge_offsets);
168
169 auto child_block_offset_it = std::prev(std::upper_bound(merge_begin, merge_end, block_offset));
170 // child_index: child node for whom `block_offset` belongs to
171 std::size_t child_index = std::distance(merge_begin, child_block_offset_it);
172 // child_block_offset: first index of the child within the block
173 auto child_block_offset = *child_block_offset_it;
174 // child_block_size: number of indices within the block partition that belong to the child
175 auto child_block_size = *std::next(child_block_offset_it) - child_block_offset;
176
177 // reconstruct index within the child ordering
178 suffix.back() = child_block_size * block_index + block_offset;
179
180 // and delegate the rest of the suffix to the child ordering
181 assert(node().degree() > child_index);
182 if constexpr (Node::isPower)
183 return node().child(child_index).size(suffix);
184 else {
185 typename Traits::SizeType _size;
186 // unfold all (tuple) children and find the child index that we want
187 Hybrid::forEach(Dune::range(node().degree()), [&](auto i){
188 if (i == child_index)
189 _size = node().child(i).size(suffix);
190 });
191 return _size;
192 }
193 }
194
195 private:
197 const Node& node() const { return static_cast<const Node&>(*this); }
198 Node& node() { return static_cast<Node&>(*this); }
199 };
200
201 } // namespace interleaved_ordering
202
203
204 template<typename DI, typename CI, typename Child, std::size_t k>
205 class PowerInterleavedOrdering
206 : public TypeTree::PowerNode<Child, k>
207 , public interleaved_ordering::Base<DI,
208 CI,
209 PowerInterleavedOrdering<DI,CI,Child,k>
210 >
211 {
212 typedef TypeTree::PowerNode<Child, k> Node;
213
214 typedef interleaved_ordering::Base<DI,
215 CI,
216 PowerInterleavedOrdering<DI,CI,Child,k>
217 > Base;
218
219 public:
220
222
230 PowerInterleavedOrdering(bool container_blocked, const InterleavedOrderingTag& ordering_tag, const typename Node::NodeStorage& children, typename Base::GFSData* gfs_data)
231 : Node(children)
232 , Base(*this,container_blocked,ordering_tag,gfs_data)
233 {}
234
235 void update()
236 {
237 for (std::size_t i = 0; i < k; ++i)
238 {
239 this->child(i).update();
240 }
241 Base::update();
242 }
243
244 std::string name() const { return "PowerInterleavedOrdering"; }
245 };
246
247
248 template<typename GFS, typename Transformation>
249 struct power_gfs_to_interleaved_ordering_descriptor
250 {
251
252 static const bool recursive = true;
253
254 template<typename TC>
255 struct result
256 {
257
258 typedef PowerInterleavedOrdering<
259 typename Transformation::DOFIndex,
260 typename Transformation::ContainerIndex,
261 TC,
262 TypeTree::StaticDegree<GFS>::value
263 > type;
264
265 typedef std::shared_ptr<type> storage_type;
266
267 };
268
269 template<typename TC>
270 static typename result<TC>::type transform(const GFS& gfs, const Transformation& t, const std::array<std::shared_ptr<TC>,TypeTree::StaticDegree<GFS>::value>& children)
271 {
272 return typename result<TC>::type(gfs.backend().blocked(gfs),gfs.orderingTag(),children,const_cast<GFS*>(&gfs));
273 }
274
275 template<typename TC>
276 static typename result<TC>::storage_type transform_storage(std::shared_ptr<const GFS> gfs, const Transformation& t, const std::array<std::shared_ptr<TC>,TypeTree::StaticDegree<GFS>::value>& children)
277 {
278 return std::make_shared<typename result<TC>::type>(gfs->backend().blocked(*gfs),gfs->orderingTag(),children,const_cast<GFS*>(gfs.get()));
279 }
280
281 };
282
283 template<typename GFS, typename Transformation>
284 power_gfs_to_interleaved_ordering_descriptor<GFS,Transformation>
285 register_power_gfs_to_ordering_descriptor(GFS*,Transformation*,InterleavedOrderingTag*);
286
287
288
289 template<typename DI, typename CI, typename... Children>
290 class CompositeInterleavedOrdering :
291 public TypeTree::CompositeNode<Children...>,
292 public interleaved_ordering::Base<DI,
293 CI,
294 CompositeInterleavedOrdering<
295 DI,
296 CI,
297 Children...
298 >
299 >
300 {
301 typedef TypeTree::CompositeNode<Children...> Node;
302
303 typedef interleaved_ordering::Base<
304 DI,
305 CI,
306 CompositeInterleavedOrdering<
307 DI,
308 CI,
309 Children...
310 >
311 > Base;
312
313 public:
315
323 CompositeInterleavedOrdering(bool backend_blocked, const InterleavedOrderingTag& ordering_tag, typename Base::GFSData* gfs_data, std::shared_ptr<Children>... children)
324 : Node(children...)
325 , Base(*this,backend_blocked,ordering_tag,gfs_data)
326 { }
327
328 std::string name() const { return "CompositeInterleavedOrdering"; }
329
330 void update()
331 {
332 TypeTree::applyToTree(*this,ordering::update_direct_children());
333 Base::update();
334 }
335 };
336
337 template<typename GFS, typename Transformation>
338 struct composite_gfs_to_interleaved_ordering_descriptor
339 {
340
341 static const bool recursive = true;
342
343 template<typename... TC>
344 struct result
345 {
346
347 typedef CompositeInterleavedOrdering<
348 typename Transformation::DOFIndex,
349 typename Transformation::ContainerIndex,
350 TC...
351 > type;
352
353 typedef std::shared_ptr<type> storage_type;
354
355 };
356
357 template<typename... TC>
358 static typename result<TC...>::type transform(const GFS& gfs, const Transformation& t, std::shared_ptr<TC>... children)
359 {
360 return typename result<TC...>::type(gfs.backend().blocked(gfs),gfs.orderingTag(),const_cast<GFS*>(&gfs),children...);
361 }
362
363 template<typename... TC>
364 static typename result<TC...>::storage_type transform_storage(std::shared_ptr<const GFS> gfs, const Transformation& t, std::shared_ptr<TC>... children)
365 {
366 return std::make_shared<typename result<TC...>::type>(gfs->backend().blocked(*gfs),gfs.orderingTag(),const_cast<GFS*>(gfs.get()),children...);
367 }
368
369 };
370
371 template<typename GFS, typename Transformation>
372 composite_gfs_to_interleaved_ordering_descriptor<GFS,Transformation>
373 register_composite_gfs_to_ordering_descriptor(GFS*,Transformation*,InterleavedOrderingTag*);
374
376 } // namespace PDELab
377} // namespace Dune
378
379#endif // DUNE_PDELAB_ORDERING_INTERLEAVEDORDERING_HH
Error related to the logical structure of an Ordering.
Definition: exceptions.hh:46
Interface for merging index spaces.
Definition: interleavedordering.hh:30
Base(Node &node, bool container_blocked, const OrderingTag &ordering_tag, typename BaseT::GFSData *gfs_data)
Construct ordering object.
Definition: interleavedordering.hh:48
std::array< std::shared_ptr< Child >, k > NodeStorage
The type used for storing the children.
Definition: powernode.hh:77
PDELab-specific exceptions.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:79
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:239
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:128
Dune namespace.
Definition: alignedallocator.hh:13
Indicate interleaved ordering of the unknowns of non-leaf grid function spaces according to a given b...
Definition: tags.hh:79
const std::vector< std::size_t > & offsets() const
Returns a list of offsets for the child blocks.
Definition: tags.hh:115
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)