Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__ntree__
34 : #define __H__UG__ntree__
35 :
36 : #include "ntree_iterator.h"
37 :
38 : namespace ug{
39 :
40 : /** The following methods have to be provided for the given vector-type:
41 : * \code
42 : * void VecSet(vector_t& vOut, real_t value); // sets all components of 'v' to 'value'.
43 : * void VecAdd(vector_t& vOut, const vector_t& v1, const vector_t& v2); // performs vOut = v1 + v2.
44 : * void VecScale(vector_t& vOut, const vector_t& v, real_t s); // performs vOut = s * v.
45 : * \endcode
46 : */
47 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
48 : struct ntree_traits
49 : {
50 : typedef int real_t;
51 : typedef int vector_t;
52 : typedef int box_t;
53 :
54 : static void calculate_center(vector_t& centerOut, const elem_t& e,
55 : const common_data_t& commonData);
56 :
57 : static void calculate_bounding_box(box_t& boxOut, const elem_t& e,
58 : const common_data_t& commonData);
59 :
60 : /// adds the given offset to box.max and subtracts it from box.min
61 : static void grow_box(box_t& boxOut, const box_t& box,
62 : const vector_t& offset);
63 :
64 : static vector_t box_diagonal(const box_t& box);
65 :
66 : static bool box_contains_point(const box_t& box, const vector_t& point);
67 :
68 : /// returns true if the given boxes intersect
69 : static bool box_box_intersection(const box_t& box1, const box_t& box2);
70 :
71 : /// returns the smallest box that contains both box1 and box2
72 : static void merge_boxes(box_t& boxOut, const box_t& box1, const box_t& box2);
73 :
74 : /// splits the given box into (2^tree_dim sub-boxes).
75 : /** The split should be performed so that the given split-point is the only common
76 : * point of all boxes.
77 : * The union of all boxes in boxesOut has to overlap the given input-box.*/
78 : static void split_box(box_t* boxesOut, const box_t& box, const vector_t& splitPoint);
79 :
80 : //// todo: the following methods should go into special traverser traits.
81 : ///** required for ContainsPointTraverser.*/
82 : static bool contains_point(const elem_t& e, const vector_t& point,
83 : const common_data_t& commonData);
84 : //
85 : ///** required for ClosestEntriesTraverser.*/
86 : // real_t distance_point_to_entry(const elem_t& e, const vector_t& point,
87 : // const elem_data_t& elemData,
88 : // const common_data_t& commonData);
89 : //
90 : ///** required for IntersectingEntriesTraverser.*/
91 : // bool ray_intersects_entry(real_t distOut, const elem_t& elem,
92 : // const vector_t& rayFrom, const vector_t& rayTo,
93 : // const elem_data_t& elemData,
94 : // const common_data_t& commonData);
95 : };
96 :
97 :
98 :
99 :
100 : struct NTreeDesc{
101 : size_t maxDepth;
102 : size_t splitThreshold;
103 : //todo: split-strategy (center-of-space, center-of-mass)
104 :
105 0 : NTreeDesc() : maxDepth(32), splitThreshold(16) {}
106 : };
107 :
108 :
109 : /// The n-tree class can be used to construct space partitioning trees of dimensions 1, 2, and 3.
110 : /** The ntree class provides spatial partitioning through binary- (tree_dim=1),
111 : * quad- (tree_dim=2), and octrees (tree_dim=3). It is designed as a template class to
112 : * support both different types of underlying math libraries and also arbitrary objects
113 : * that shall be partitioned. The element type is specified through TElem. It
114 : * does not have to fulfill any concepts but should be lightweight, since it is
115 : * copied during tree creation.
116 : *
117 : * The tree furthermore features a common_data object whose type is determined
118 : * through a template argument again. This common_data object is passed to
119 : * traversers and traits. Required objects such as data-accessors can be stored
120 : * in this common_data object as required by a specialization of ntree_traits.
121 : * It does not have to fulfill any concepts in regard to the ntree itself.
122 : *
123 : * Usage:
124 : * - Add elements to your tree-instance through the 'add_element' method.
125 : * - Call 'rebalance' once all elements have been added.
126 : * - Use traversers to query the tree for elements, e.g.
127 : * 'Traverser_FindContainingElement', 'Traverser_FindElementsInIntersectingNodes',
128 : * or 'Traverser_RayElementIntersection'.
129 : *
130 : * \note Whether a node is refined into child-nodes during 'rebalance' is
131 : * determined by the 'splitThreshold' of the associated 'NTreeDesc'
132 : * instance. It can be adjusted through 'set_desc'.
133 : * Only if a node contains more than 'splitThreshold' elements, a split
134 : * of that node will be performed. Default is 16.
135 : *
136 : * \note Each tree has a maximum tree depth. It defaults to 32 and can be
137 : * adjusted through the 'set_desc' method. If this tree depth is
138 : * reached, the tree won't be splitted any further and a warning
139 : * is printed. Note that a tree of this depth often results from a bad
140 : * underlying geometry and may lead to performance issues.
141 : *
142 : * \param tree_dim Dimension of the tree: binary-trees (tree_dim=1),
143 : * quad-trees (tree_dim=2), and octrees (tree_dim=3)
144 : * are supported.
145 : *
146 : * \param world_dim Dimension of the space in which the tree is embedded.
147 : * This primarily affects underlying mathematical structures
148 : * (i.e. vectors). Please note that world_dim has to be at
149 : * least as large as tree_dim.
150 : *
151 : * \param TElem The element type for which the tree is constructed. It
152 : * should be lightweight, since it is copied during tree
153 : * creation. There are no special concepts that TElem has
154 : * to fullfill.
155 : *
156 : * \param TCommonData User provided data that is stored in the tree (one instance only)
157 : * but not used by the tree. It is passed to functions in
158 : * ntree_traits. The concrete type for TCommonData depends on
159 : * the concrete ntree_traits that shall be used.
160 : */
161 : template <int tree_dim, int world_dim, class TElem, class TCommonData>
162 0 : class ntree
163 : {
164 : private:
165 : struct Entry;
166 :
167 : public:
168 : typedef TElem elem_t;
169 : typedef TCommonData common_data_t;
170 : typedef ntree_traits<tree_dim, world_dim, elem_t, common_data_t> traits;
171 : typedef typename traits::real_t real_t;
172 : typedef typename traits::vector_t vector_t;
173 : typedef typename traits::box_t box_t;
174 : typedef const_ntree_element_iterator<elem_t, Entry> elem_iterator_t;
175 :
176 : ntree();
177 :
178 : void clear();
179 :
180 : /// enabled or disable warning messages
181 : /** warnings are enabled by default. If a problem is detected and warnings
182 : * are enabled, a warning message will be written to stdout.*/
183 : void enable_warnings(bool enable) {m_warningsEnabled = enable;}
184 : bool warnings_enabled () const {return m_warningsEnabled;}
185 :
186 : /// sets the balancing-parameters of the tree.
187 : /** \note The method has no effect until the next call to 'rebalance',
188 : * i.e., it will not affect an existing tree if one has already
189 : * been built.*/
190 : void set_desc(const NTreeDesc& desc);
191 :
192 : /// returns the balancing-parameters of the tree.
193 : const NTreeDesc& desc() const;
194 :
195 : /// sets the common-data which the tree passes on to callback methods
196 : void set_common_data(const common_data_t& commonData);
197 :
198 : /// returns the common-data stored in the tree
199 : const common_data_t& common_data() const;
200 :
201 : /// returns true if the tree is empty
202 : bool empty() const;
203 :
204 : /// returns the number of entries in the tree (delayed entries excluded)
205 : size_t size() const;
206 :
207 : /// returns the number of elements which have been added but are not yet accessible in the tree.
208 : /** \note delayed elements are inserted on a call to rebalance.*/
209 : size_t num_delayed_elements() const;
210 :
211 : /// adds an element to the tree.
212 : /** The element will only be scheduled for insertion but won't be inserted
213 : * until rebalance is called.
214 : * \todo Allow optional on-the-fly insertion of elements.*/
215 : void add_element(const elem_t& elem);
216 :
217 : /// rebalances the whole tree
218 : /** The method returns false if some error occurred during rebalancing.*/
219 : void rebalance();
220 :
221 : /// returns the total number of nodes in the tree
222 : size_t num_nodes() const;
223 :
224 : /// returns the number of children of a node
225 : size_t num_child_nodes(size_t nodeId) const;
226 :
227 : /// returns an array of child-id's for the given node
228 : /** Use ntree::num_children on the node to retrieve the length of the returned array.*/
229 : const size_t* child_node_ids(size_t nodeId) const;
230 :
231 : /// returns an iterator to the first element of a given node
232 : elem_iterator_t elems_begin(size_t nodeId) const;
233 :
234 : /// returns an iterator to the end of the element-sequence of a given node
235 : /** this iterator points one element behind the last element of the sequence.*/
236 : elem_iterator_t elems_end(size_t nodeId) const;
237 :
238 : /// returns the number of elements that the given node contains
239 : size_t num_elements(size_t nodeId) const;
240 :
241 : /// returns the number tree-level in which the node is located
242 : size_t level(size_t nodeId) const;
243 :
244 : /// returns the smallest box which contains all elements of the given node
245 : const box_t& bounding_box(size_t nodeId) const;
246 :
247 : private:
248 :
249 : /// clear only the nodes
250 : void clear_nodes ();
251 :
252 : /// static template implementation to raise n to the power exponent
253 : template <size_t n, size_t exponent>
254 : struct pow;
255 :
256 : template <size_t n>
257 : struct pow<n, 0> {static const size_t val = 1;};
258 :
259 : template <size_t n, size_t exponent>
260 : struct pow {static const size_t val = n * pow<n, exponent - 1>::val;};
261 :
262 : /// the number of children each non-leaf node has
263 : static const size_t s_numChildren = pow<2, tree_dim>::val;
264 :
265 : /// marks an index as invalid
266 : static const size_t s_invalidIndex = -1;
267 :
268 :
269 : /// An Entry combines an element with the index to the next entry in a leaf-node's entry list.
270 : /** Note that exactly one 'Entry' object per element exists. Since 'Entry'
271 : * also serves as a linked list, this means that an element can only be contained
272 : * in one node at a time.*/
273 : struct Entry{
274 : elem_t elem;
275 : /** index into m_entries. s_invalidIndex: no next entry.
276 : * Used to create a linked list for each leaf-node.*/
277 : size_t nextEntryInd;
278 :
279 0 : Entry(const elem_t& e) :
280 0 : elem(e), nextEntryInd(s_invalidIndex) {}
281 : };
282 :
283 :
284 : /** The tree is built as a hierarchy of nodes.
285 : * Leaf nodes contain entries.*/
286 : struct Node{
287 : size_t childNodeInd[s_numChildren]; /// < index into m_nodes. s_invalidIndex: no child node.
288 : size_t firstEntryInd; ///< index into m_entries. s_invalidIndex: no entry
289 : size_t lastEntryInd; ///< index into m_entries. s_invalidIndex: no entry
290 : size_t numEntries; ///< number of entries in the node
291 : size_t level;
292 : box_t tightBox; ///< tight bounding box - disjunct partition of the root box
293 : box_t looseBox; ///< loose bounding box - contains all bounding boxes of its entries
294 :
295 0 : Node() : firstEntryInd(s_invalidIndex), lastEntryInd(s_invalidIndex),
296 0 : numEntries(0), level(s_invalidIndex)
297 : {
298 : tightBox = looseBox;
299 0 : for(size_t i = 0; i < s_numChildren; ++i)
300 0 : childNodeInd[i] = s_invalidIndex;
301 : }
302 :
303 0 : Node(const Node& n) : firstEntryInd(n.firstEntryInd), lastEntryInd(n.lastEntryInd),
304 0 : numEntries(n.numEntries), level(n.level),
305 : tightBox(n.tightBox), looseBox(n.looseBox)
306 : {
307 0 : for(size_t i = 0; i < s_numChildren; ++i)
308 0 : childNodeInd[i] = n.childNodeInd[i];
309 0 : }
310 : };
311 :
312 : /// splits a node into 2^tree_dim child nodes and assigns entries to those children.
313 : /** If the node-threshold of a child node is surpassed, then the child will
314 : * be splitted recursively.
315 : * Make sure that the given node is a leaf node and thus hasn't got children.*/
316 : void split_leaf_node(size_t nodeIndex);
317 :
318 : /// returns an index to the leaf-node which contains the given point
319 : /** returns s_invalidIndex if no matching node was found.
320 : * Checks the point against the tight bounding-box of each node.*/
321 : size_t find_leaf_node(const vector_t& point, size_t curNode = 0);
322 :
323 : /// adds an entry to the given node
324 : void add_entry_to_node(Node& node, size_t entryInd);
325 :
326 : /// updates the loose bounding box of the given node
327 : void update_loose_bounding_box(Node& node);
328 :
329 : /// calculates the center of mass of a given node
330 : vector_t calculate_center_of_mass(Node& node);
331 :
332 : NTreeDesc m_desc;
333 : common_data_t m_commonData;
334 : std::vector<Node> m_nodes; ///< m_nodes[0] is always considered to be the root node.
335 : std::vector<Entry> m_entries;
336 : size_t m_numDelayedElements;
337 : bool m_warningsEnabled;
338 : };
339 :
340 :
341 : }// end of namespace
342 :
343 :
344 : ////////////////////////////////////////
345 : // include implementation
346 : #include "ntree_impl.hpp"
347 :
348 : #endif
|