1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_ORDERING_UTILITY_HH
5 #define DUNE_PDELAB_ORDERING_UTILITY_HH
6 
7 #include <vector>
8 #include <bitset>
9 
10 #include <dune/pdelab/common/dofindex.hh>
11 #include <dune/pdelab/common/globaldofindex.hh>
12 #include <dune/pdelab/ordering/transformations.hh>
13 
14 namespace Dune {
15   namespace PDELab {
16 
17     //! \addtogroup Ordering
18     //! \{
19 
20     //! Index merging algorithm for global orderings.
21     struct MergeMode
22     {
23 
24       enum type {
25         lexicographic, //!< Lexicographically ordered ([i1,i2],[j1,j2] -> [i1,i2,j1,j2]).
26         interleaved  //!< Indices are interleaved according to a user-supplied pattern ([i1,i2],[j1,j2] -> [i1,j1,i2,j2]).
27       };
28 
29     };
30 
31     //! Information about order semantics on multi-indices
32     enum class MultiIndexOrder {
33       //! indices are ordered from inner to outer container: {inner,...,outer}
34       Inner2Outer,
35       //! indices are ordered from outer to inner container: {outer,...,inner}
36       Outer2Inner,
37     };
38 
39 
40 #ifndef DOXYGEN
41 
42     namespace ordering {
43 
44       // This is an implementation detail of the composite orderings, no need to confuse our users!
45       struct update_direct_children
46         : public TypeTree::DirectChildrenVisitor
47         , public TypeTree::DynamicTraversal
48       {
49 
50         template<typename GFS, typename Child, typename TreePath, typename ChildIndex>
afterChildDune::PDELab::ordering::update_direct_children51         void afterChild(const GFS& gfs, Child& child, TreePath tp, ChildIndex childIndex) const
52         {
53           child.update();
54         }
55 
56       };
57 
58     } // end namespace ordering
59 
60 #endif // DOXYGEN
61 
62 
63     struct DefaultDOFIndexAccessor
64     {
65 
66       template<typename DOFIndex, typename SizeType, typename IndexType>
67       static typename std::enable_if<
68         std::is_integral<IndexType>::value
69         >::type
storeDune::PDELab::DefaultDOFIndexAccessor70       store(DOFIndex& dof_index, const GeometryType& gt, SizeType entity_index, IndexType tree_index)
71       {
72         dof_index.clear();
73         dof_index.entityIndex()[0] = GlobalGeometryTypeIndex::index(gt);
74         dof_index.entityIndex()[1] = entity_index;
75         dof_index.treeIndex().push_back(tree_index);
76       }
77 
78       template<typename DOFIndex, typename SizeType, typename IndexType>
79       static typename std::enable_if<
80         !std::is_integral<IndexType>::value
81         >::type
storeDune::PDELab::DefaultDOFIndexAccessor82       store(DOFIndex& dof_index, const GeometryType& gt, SizeType entity_index, IndexType tree_index)
83       {
84         dof_index.entityIndex()[0] = GlobalGeometryTypeIndex::index(gt);
85         dof_index.entityIndex()[1] = entity_index;
86         dof_index.treeIndex() = tree_index;
87       }
88 
89       template<typename DOFIndex, typename SizeType, typename IndexType>
90       static typename std::enable_if<
91         std::is_integral<IndexType>::value
92         >::type
storeDune::PDELab::DefaultDOFIndexAccessor93       store(DOFIndex& dof_index, SizeType gt_index, SizeType entity_index, IndexType tree_index)
94       {
95         dof_index.clear();
96         dof_index.entityIndex()[0] = gt_index;
97         dof_index.entityIndex()[1] = entity_index;
98         dof_index.treeIndex().push_back(tree_index);
99       }
100 
101       template<typename DOFIndex, typename SizeType, typename IndexType>
102       static typename std::enable_if<
103         !std::is_integral<IndexType>::value
104         >::type
storeDune::PDELab::DefaultDOFIndexAccessor105       store(DOFIndex& dof_index, SizeType gt_index, SizeType entity_index, IndexType tree_index)
106       {
107         dof_index.entityIndex()[0] = gt_index;
108         dof_index.entityIndex()[1] = entity_index;
109         dof_index.treeIndex() = tree_index;
110       }
111 
112 
113       struct GeometryIndex
114       {
115 
116         template<typename Index>
geometryTypeDune::PDELab::DefaultDOFIndexAccessor::GeometryIndex117         static std::size_t geometryType(const Index& geometry_index)
118         {
119           return geometry_index[0];
120         }
121 
122         template<typename Index>
entityIndexDune::PDELab::DefaultDOFIndexAccessor::GeometryIndex123         static std::size_t entityIndex(const Index& geometry_index)
124         {
125           return geometry_index[1];
126         }
127 
128         template<typename Index, typename SizeType>
storeDune::PDELab::DefaultDOFIndexAccessor::GeometryIndex129         static void store(Index& index, const GeometryType& gt, SizeType entity_index)
130         {
131           index[0] = GlobalGeometryTypeIndex::index(gt);
132           index[1] = entity_index;
133         }
134 
135         template<typename Index, typename SizeType>
storeDune::PDELab::DefaultDOFIndexAccessor::GeometryIndex136         static void store(Index& index, SizeType geometry_index, SizeType entity_index)
137         {
138           index[0] = geometry_index;
139           index[1] = entity_index;
140         }
141 
142       };
143 
144       template<typename DOFIndex>
geometryTypeDune::PDELab::DefaultDOFIndexAccessor145       static std::size_t geometryType(const DOFIndex& dof_index)
146       {
147         return GeometryIndex::geometryType(dof_index.entityIndex());
148       }
149 
150       template<typename DOFIndex>
entityIndexDune::PDELab::DefaultDOFIndexAccessor151       static std::size_t entityIndex(const DOFIndex& dof_index)
152       {
153         return GeometryIndex::entityIndex(dof_index.entityIndex());
154       }
155 
156     };
157 
158     struct SimpleDOFIndexAccessor
159     {
160 
161       template<typename DOFIndex, typename SizeType>
storeDune::PDELab::SimpleDOFIndexAccessor162       static void store(DOFIndex& dof_index, const GeometryType& gt, SizeType entity_index, SizeType tree_index)
163       {
164         dof_index = entity_index;
165       }
166 
167     };
168 
169 
170     template<typename DI, typename CI, MultiIndexOrder CIOrder = MultiIndexOrder::Inner2Outer>
171     struct SimpleOrderingTraits
172     {
173 
174       typedef DI DOFIndex;
175 
176       typedef CI ContainerIndex;
177 
178       typedef std::size_t SizeType;
179 
180       typedef DefaultDOFIndexAccessor DOFIndexAccessor;
181 
182       //! Inform about ContainerIndex multi-index order semantics
183       static constexpr MultiIndexOrder ContainerIndexOrder = CIOrder;
184     };
185 
186 
187     template<typename SizeType_, typename CI, MultiIndexOrder CIOrder>
188     struct SimpleOrderingTraits<SimpleDOFIndex<SizeType_>,CI,CIOrder>
189     {
190 
191       typedef SimpleDOFIndex<SizeType_> DOFIndex;
192 
193       typedef CI ContainerIndex;
194 
195       typedef SizeType_ SizeType;
196 
197       typedef SimpleDOFIndexAccessor DOFIndexAccessor;
198 
199       //! Inform about ContainerIndex multi-index order semantics
200       static constexpr MultiIndexOrder ContainerIndexOrder = CIOrder;
201     };
202 
203     template <typename DI, typename CI,
204               MultiIndexOrder CIOrder = MultiIndexOrder::Inner2Outer>
205     struct OrderingTraits : public SimpleOrderingTraits<DI, CI, CIOrder> {
206 
207       // The maximum dimension supported (length of bitsets)
208       // 32 dimensions should probably be fine for now... ;-)
209       static const std::size_t max_dim = 32;
210 
211       typedef std::bitset<max_dim> CodimFlag;
212 
213       typedef typename DI::TreeIndex TreeIndex;
214 
215       typedef typename DI::View DOFIndexView;
216       typedef typename DI::View::TreeIndex TreeIndexView;
217 
218       typedef typename DI::size_type SizeType;
219       typedef typename DI::size_type size_type;
220     };
221 
222     template <typename ES, typename DI, typename CI,
223               MultiIndexOrder CIOrder = MultiIndexOrder::Outer2Inner>
224     struct LocalOrderingTraits : public OrderingTraits<DI, CI, CIOrder> {
225 
226       using EntitySet = ES;
227       using GridView = typename ES::GridView;
228     };
229 
230     template<typename ES, typename DI, typename CI>
231     struct GridViewOrderingTraits
232       : public LocalOrderingTraits<ES,DI,CI>
233     {
234 
235       typedef typename DI::EntityIndex EntityIndex;
236       typedef typename DI::View::EntityIndex EntityIndexView;
237 
238     };
239 
240 
241     template<typename DI, typename CI>
242     class VirtualOrderingBase
243     {
244     public:
245 
246       typedef OrderingTraits<DI,CI> Traits;
247 
VirtualOrderingBase()248       VirtualOrderingBase() {}
249       virtual ~VirtualOrderingBase() = default;
250 
251       virtual void map_index_dynamic(typename Traits::DOFIndexView di, typename Traits::ContainerIndex& ci) const = 0;
252     };
253 
254 
255     template<typename child_type>
256     struct extract_child_bases
257       : public TypeTree::DirectChildrenVisitor
258       , public TypeTree::DynamicTraversal
259     {
260 
261       template<typename Node, typename Child, typename TreePath, typename ChildIndex>
afterChildDune::PDELab::extract_child_bases262       void afterChild(const Node& node, Child& child, TreePath tp, ChildIndex child_index)
263       {
264         extract_child(node,child,child_index);
265       }
266 
267       template<typename Node, typename Child, typename ChildIndex>
268       typename std::enable_if<Node::has_dynamic_ordering_children>::type
extract_childDune::PDELab::extract_child_bases269       extract_child(const Node& node, Child& child, ChildIndex child_index)
270       {
271         _children[child_index] = &child;
272       }
273 
274       template<typename Node, typename Child, typename ChildIndex>
275       typename std::enable_if<!Node::has_dynamic_ordering_children>::type
extract_childDune::PDELab::extract_child_bases276       extract_child(const Node& node, Child& child, ChildIndex child_index)
277       {
278       }
279 
extract_child_basesDune::PDELab::extract_child_bases280       extract_child_bases(std::vector<child_type*>& children)
281         : _children(children)
282       {}
283 
284     private:
285       std::vector<child_type*>& _children;
286 
287     };
288 
289 
290     //! Dummy iterator type over DOF indices.
291     /**
292      * This dummy iterator is used to support omitting the calculation
293      * of DOFIndex values in the per-entity index lookup methods of
294      * orderings. By defining all operations performed on the DOFIndex
295      * iterator and its value by this methods as no-ops, we can reuse the
296      * combined implementation mapping both DOFIndex and ContainerIndex for
297      * the (much more common) case of only having to map the ContainerIndex
298      * values.
299      */
300     struct DummyDOFIndexIterator
301     {
302 
303       typedef std::size_t size_type;
304 
operator ++Dune::PDELab::DummyDOFIndexIterator305       DummyDOFIndexIterator& operator++()
306       {
307         return *this;
308       }
309 
operator +=Dune::PDELab::DummyDOFIndexIterator310       DummyDOFIndexIterator& operator+=(size_type i)
311       {
312         return *this;
313       }
314 
operator *Dune::PDELab::DummyDOFIndexIterator315       DummyDOFIndexIterator& operator*()
316       {
317         return *this;
318       }
319 
operator ->Dune::PDELab::DummyDOFIndexIterator320       DummyDOFIndexIterator* operator->()
321       {
322         return this;
323       }
324 
treeIndexDune::PDELab::DummyDOFIndexIterator325       DummyDOFIndexIterator& treeIndex()
326       {
327         return *this;
328       }
329 
operator ==Dune::PDELab::DummyDOFIndexIterator330       bool operator==(const DummyDOFIndexIterator& r) const
331       {
332         return true;
333       }
334 
operator !=Dune::PDELab::DummyDOFIndexIterator335       bool operator!=(const DummyDOFIndexIterator& r) const
336       {
337         return !operator==(r);
338       }
339 
push_backDune::PDELab::DummyDOFIndexIterator340       void push_back(size_type i)
341       {}
342 
343     };
344 
345     /**
346      * @brief Adapter to create a size provider from an ordering
347      * @details This adapter is meant to be used in allocation and
348      *   resizing of vectors containers.
349      *   In particular, this adapter is needed because the ordering library give
350      *   sizes for multi-indices ordered with Inner2Outer semantis, while
351      *   resizing algorithms are faster and easier when using Outer2Inner
352      *   semantics.
353      *
354      *   - This class makes type erasure on the ordering type.
355      *   - This class has value semantics.
356      *   - This class always receives container indices with Outer2Inner order
357      *
358      * @tparam SizeType_ return type of the size method
359      * @tparam ContainerIndex_ argument type of the size method
360      * @tparam OriginOrder enum with MultiIndexOrder semantics of the origin ordering
361      */
362     template <class Size, class ContainerIndex_, MultiIndexOrder OriginOrder>
363     struct SizeProviderAdapter {
364 
365       /**
366        * @brief Construct a new Size Provider Adapter object
367        *
368        * @tparam Ordering  The type of the ordering to adapt
369        * @param ordering   A shared pointer to the ordering
370        */
371       template <class Ordering>
SizeProviderAdapterDune::PDELab::SizeProviderAdapter372       SizeProviderAdapter(const std::shared_ptr<const Ordering> &ordering)
373           : _size_provider([=](const ContainerIndex_ &partial_multiindex) {
374               return ordering->size(partial_multiindex);
375             }) {
376         static_assert(Ordering::Traits::ContainerIndexOrder == OriginOrder);
377       }
378 
379       //! Partial MultiIndex of a ContainerIndex
380       using ContainerIndex = ContainerIndex_;
381 
382       //! Partial MultiIndex of a ContainerIndex
383       using SizePrefix = ContainerIndex;
384 
385       //! Type that refers to the size of containers
386       using SizeType = Size;
387 
388       //! Inform about ContainerIndex multi-index order semantics
389       static constexpr MultiIndexOrder ContainerIndexOrder = MultiIndexOrder::Outer2Inner;
390 
391       /**
392        * @brief Gives the size for a given prefix
393        * @param prefix  MultiIndex with a partial path to a container
394        * @return Traits::SizeType  The size required for such a path
395        */
sizeDune::PDELab::SizeProviderAdapter396       SizeType size(const SizePrefix &prefix) const {
397         if constexpr (OriginOrder == MultiIndexOrder::Inner2Outer) {
398           // reversing Outer2Inner prefix into a Inner2Outer suffix
399           ContainerIndex suffix;
400           suffix.resize(prefix.size());
401           std::reverse_copy(prefix.begin(), prefix.end(),suffix.begin());
402           // forward size request to ordering with new Inner2Outer suffix
403           return _size_provider(suffix);
404         } else {
405           // prefix is already Outer2Inner, forward to size provider
406           return _size_provider(prefix);
407         }
408       }
409 
410     private:
411 
412       const std::function<SizeType(const ContainerIndex &)> _size_provider;
413     };
414 
415     //! template deduction guide for orderings
416     template<class Ordering>
417     SizeProviderAdapter(const std::shared_ptr<const Ordering>& ordering)
418         -> SizeProviderAdapter<typename Ordering::Traits::SizeType,
419                               typename Ordering::Traits::ContainerIndex,
420                               Ordering::Traits::ContainerIndexOrder>;
421 
422 
423    //! \} group Ordering
424   } // namespace PDELab
425 } // namespace Dune
426 
427 #endif // DUNE_PDELAB_ORDERING_UTILITY_HH
428