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_impl__
34 : #define __H__UG__ntree_impl__
35 :
36 : #include <cassert>
37 : #include "ntree.h"
38 :
39 : namespace ug{
40 :
41 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
42 0 : ntree<tree_dim, world_dim, elem_t, common_data_t>::
43 : ntree() :
44 0 : m_warningsEnabled (true)
45 : {
46 0 : m_nodes.resize(1);
47 0 : }
48 :
49 :
50 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
51 : void ntree<tree_dim, world_dim, elem_t, common_data_t>::
52 : clear_nodes()
53 : {
54 : m_nodes.clear();
55 0 : m_nodes.resize(1);
56 0 : m_nodes[0].level = 0;
57 : }
58 :
59 :
60 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
61 0 : void ntree<tree_dim, world_dim, elem_t, common_data_t>::
62 : clear()
63 : {
64 : clear_nodes();
65 : m_entries.clear();
66 0 : m_numDelayedElements = 0;
67 0 : }
68 :
69 :
70 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
71 : void ntree<tree_dim, world_dim, elem_t, common_data_t>::
72 : set_desc(const NTreeDesc& desc)
73 : {
74 : m_desc = desc;
75 : }
76 :
77 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
78 : const NTreeDesc& ntree<tree_dim, world_dim, elem_t, common_data_t>::
79 : desc() const
80 : {
81 : return m_desc;
82 : }
83 :
84 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
85 : void ntree<tree_dim, world_dim, elem_t, common_data_t>::
86 : set_common_data(const common_data_t& commonData)
87 : {
88 0 : m_commonData = commonData;
89 : }
90 :
91 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
92 : const common_data_t& ntree<tree_dim, world_dim, elem_t, common_data_t>::
93 : common_data() const
94 : {
95 0 : return m_commonData;
96 : }
97 :
98 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
99 : bool ntree<tree_dim, world_dim, elem_t, common_data_t>::
100 : empty() const
101 : {
102 : return size() == 0;
103 : }
104 :
105 :
106 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
107 : size_t ntree<tree_dim, world_dim, elem_t, common_data_t>::
108 : size() const
109 : {
110 : return m_entries.size() - m_numDelayedElements;
111 : }
112 :
113 :
114 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
115 : size_t ntree<tree_dim, world_dim, elem_t, common_data_t>::
116 : num_delayed_elements() const
117 : {
118 : return m_numDelayedElements;
119 : }
120 :
121 :
122 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
123 : void ntree<tree_dim, world_dim, elem_t, common_data_t>::
124 : add_element(const elem_t& elem)
125 : {
126 : // size_t entryInd = m_entries.size();
127 0 : m_entries.push_back(Entry(elem));
128 :
129 0 : ++m_numDelayedElements;
130 :
131 : //todo If we allowed for immediate on-the-fly insertion of elements, the following
132 : // code could serve as an implementation base. Note, that it is not yet complete.
133 : // One would have to update parent boxes if a child box changes etc.
134 : // It would make sense to also allow for the specification of a
135 : // minimum root-bounding-box, since else most insertions would trigger a rebalance.
136 : // {
137 : // vector_t center;
138 : // traits::calculate_center(center, elem, m_commonData);
139 : // size_t nodeInd = find_leaf_node(center);
140 : // if(nodeInd == s_invalidIndex){
141 : // ++m_numDelayedElements;
142 : // rebalance();
143 : // }
144 : // else{
145 : // Node& node = m_nodes[nodeInd];
146 : // add_entry_to_node(node, entryInd);
147 : // if(node.numEntries >= m_desc.splitThreshold)
148 : // split_leaf_node(nodeInd);
149 : // else{
150 : // //todo the update could be much more efficient (simply unite the
151 : // // old bounding box with the new element's bounding box.
152 : // update_loose_bounding_box(node);
153 : // }
154 : // }
155 : // }
156 : }
157 :
158 :
159 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
160 0 : void ntree<tree_dim, world_dim, elem_t, common_data_t>::
161 : rebalance()
162 : {
163 : // push all elements into the root node, calculate its bounding box
164 : // and call split_leaf_node if the element threshold is surpassed.
165 : clear_nodes();
166 : Node& root = m_nodes.back();
167 0 : if(!m_entries.empty()){
168 0 : root.firstEntryInd = 0;
169 0 : root.lastEntryInd = m_entries.size() - 1;
170 0 : root.numEntries = m_entries.size();
171 0 : for(size_t i = 0; i < m_entries.size(); ++i)
172 0 : m_entries[i].nextEntryInd = i+1;
173 0 : m_entries.back().nextEntryInd = s_invalidIndex;
174 :
175 : // the tight bounding box and the loose bounding box of the root node
176 : // should be the same.
177 0 : update_loose_bounding_box(root);
178 : root.tightBox = root.looseBox;
179 :
180 0 : if(root.numEntries >= m_desc.splitThreshold)
181 0 : split_leaf_node(0);
182 : }
183 0 : }
184 :
185 :
186 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
187 0 : void ntree<tree_dim, world_dim, elem_t, common_data_t>::
188 : split_leaf_node(size_t nodeIndex)
189 : {
190 :
191 0 : if(m_nodes[nodeIndex].numEntries <= 1)
192 0 : return;
193 :
194 0 : if(m_nodes[nodeIndex].level >= m_desc.maxDepth){
195 0 : if(m_warningsEnabled){
196 0 : UG_LOG("WARNING in ntree::split_leaf_node(): maximum tree depth "
197 : << m_desc.maxDepth << " reached. No further splits are performed for "
198 : " this node. Note that too many elements per node may lead to performance issues.\n"
199 : << " Number of elements in this node: " << m_nodes[nodeIndex].numEntries << std::endl
200 : << " Corner coordinates of this node: " << m_nodes[nodeIndex].tightBox << std::endl);
201 : }
202 0 : return;
203 : }
204 :
205 0 : if(m_nodes[nodeIndex].childNodeInd[0] != s_invalidIndex)
206 : return;
207 :
208 : const size_t firstChild = m_nodes.size();
209 0 : m_nodes.resize(firstChild + s_numChildren);
210 :
211 : // ATTENTION: Be careful not to resize m_nodes while using node, since this would invalidate the reference!
212 : Node& node = m_nodes[nodeIndex];
213 :
214 : // calculate center of mass and use the traits class to split the box of
215 : // the current node into 's_numChildren' child boxes. Each child box thereby
216 : // spanned by one of the corners of the original box and 'centerOfMass'.
217 0 : vector_t centerOfMass = calculate_center_of_mass(node);
218 0 : box_t childBoxes[s_numChildren];
219 0 : traits::split_box(childBoxes, node.tightBox, centerOfMass);
220 :
221 : // iterate over all entries in the current node and assign them to child nodes.
222 : size_t numEntriesAssigned = 0;
223 0 : for(size_t entryInd = node.firstEntryInd; entryInd != s_invalidIndex;){
224 : Entry& entry = m_entries[entryInd];
225 0 : size_t nextEntryInd = entry.nextEntryInd;
226 :
227 : size_t i_child;
228 : vector_t center;
229 0 : traits::calculate_center(center, entry.elem, m_commonData);
230 0 : for(i_child = 0; i_child < s_numChildren; ++i_child){
231 0 : if(traits::box_contains_point(childBoxes[i_child], center)){
232 0 : add_entry_to_node(m_nodes[firstChild + i_child], entryInd);
233 0 : ++numEntriesAssigned;
234 0 : break;
235 : }
236 : }
237 : /*-- For debugging only: --*
238 : if(i_child == s_numChildren){
239 : UG_LOG ("ERROR in ntree::split_leaf_node(): Element with center @ " << center
240 : << " does not belong to any child of the box " << node.tightBox << std::endl);
241 : }
242 : *--*/
243 :
244 : entryInd = nextEntryInd;
245 : }
246 :
247 : // all elements of the current box now should be assigned to child boxes.
248 : // we thus clear element lists and entry-count from the current node.
249 0 : UG_COND_THROW(numEntriesAssigned != node.numEntries, "Couldn't find a matching "
250 : "child node for some elements during split_leaf_node in "
251 : "ntree::split_leaf_node");
252 :
253 0 : node.firstEntryInd = node.lastEntryInd = s_invalidIndex;
254 0 : node.numEntries = 0;
255 :
256 0 : for(size_t i_child = 0; i_child < s_numChildren; ++i_child){
257 0 : node.childNodeInd[i_child] = firstChild + i_child;
258 : Node& childNode = m_nodes[firstChild + i_child];
259 0 : childNode.level = node.level + 1;
260 0 : childNode.tightBox = childBoxes[i_child];
261 0 : update_loose_bounding_box(childNode);
262 : }
263 :
264 : // since split_leaf_node resizes m_nodes and since this invalidates any references
265 : // to m_nodes, we perform the recursion in a last step
266 0 : for(size_t i_child = 0; i_child < s_numChildren; ++i_child){
267 0 : size_t childNodeInd = firstChild + i_child;
268 0 : if(m_nodes[childNodeInd].numEntries >= m_desc.splitThreshold)
269 0 : split_leaf_node(childNodeInd);
270 : }
271 : }
272 :
273 :
274 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
275 : size_t ntree<tree_dim, world_dim, elem_t, common_data_t>::
276 : find_leaf_node(const vector_t& point, size_t curNode)
277 : {
278 : Node& n = m_nodes[curNode];
279 : if(traits::box_contains_point(n.tightBox, point)){
280 : if(n.childNodeInd[0] != s_invalidIndex){
281 : for(size_t i = 0; i < s_numChildren; ++i){
282 : size_t result = find_leaf_node(point, n.childNodeInd[i]);
283 : if(result != s_invalidIndex)
284 : return result;
285 : }
286 : }
287 : else{
288 : return curNode;
289 : }
290 : }
291 : return s_invalidIndex;
292 : }
293 :
294 :
295 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
296 : void ntree<tree_dim, world_dim, elem_t, common_data_t>::
297 : add_entry_to_node(Node& node, size_t entryInd)
298 : {
299 0 : if(node.firstEntryInd == s_invalidIndex){
300 : assert(node.numEntries == 0);
301 0 : node.firstEntryInd = node.lastEntryInd = entryInd;
302 : }
303 : else{
304 0 : m_entries[node.lastEntryInd].nextEntryInd = entryInd;
305 0 : node.lastEntryInd = entryInd;
306 : }
307 :
308 0 : m_entries[entryInd].nextEntryInd = s_invalidIndex;
309 0 : ++node.numEntries;
310 : }
311 :
312 :
313 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
314 0 : void ntree<tree_dim, world_dim, elem_t, common_data_t>::
315 : update_loose_bounding_box(Node& node)
316 : {
317 0 : size_t entryInd = node.firstEntryInd;
318 0 : if(entryInd == s_invalidIndex){
319 : node.looseBox = node.tightBox;
320 0 : return;
321 : }
322 :
323 : Entry& entry = m_entries[entryInd];
324 0 : traits::calculate_bounding_box(node.looseBox, entry.elem, m_commonData);
325 :
326 0 : entryInd = entry.nextEntryInd;
327 0 : while(entryInd != s_invalidIndex){
328 : Entry& entry = m_entries[entryInd];
329 : box_t nbox;
330 0 : traits::calculate_bounding_box(nbox, entry.elem, m_commonData);
331 0 : traits::merge_boxes(node.looseBox, node.looseBox, nbox);
332 0 : entryInd = entry.nextEntryInd;
333 : }
334 :
335 : //todo: Solve box-growing in a more portable way!
336 : typename traits::vector_t offset;
337 : offset = SMALL;
338 : traits::grow_box(node.looseBox, node.looseBox, offset);
339 : }
340 :
341 :
342 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
343 : typename ntree<tree_dim, world_dim, elem_t, common_data_t>::vector_t
344 0 : ntree<tree_dim, world_dim, elem_t, common_data_t>::
345 : calculate_center_of_mass(Node& node)
346 : {
347 : vector_t centerOfMass;
348 : VecSet(centerOfMass, 0);
349 :
350 0 : for(size_t entryInd = node.firstEntryInd; entryInd != s_invalidIndex;){
351 : Entry& entry = m_entries[entryInd];
352 :
353 : vector_t center;
354 0 : traits::calculate_center(center, entry.elem, m_commonData);
355 : VecAdd(centerOfMass, centerOfMass, center);
356 0 : entryInd = entry.nextEntryInd;
357 : }
358 :
359 0 : if(node.numEntries > 0)
360 0 : VecScale(centerOfMass, centerOfMass, (real_t)1 / (real_t)node.numEntries);
361 :
362 0 : return centerOfMass;
363 : }
364 :
365 :
366 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
367 : size_t ntree<tree_dim, world_dim, elem_t, common_data_t>::
368 : num_nodes() const
369 : {
370 : return m_nodes.size();
371 : }
372 :
373 :
374 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
375 : size_t ntree<tree_dim, world_dim, elem_t, common_data_t>::
376 : num_child_nodes(size_t nodeId) const
377 : {
378 : assert(nodeId < m_nodes.size());
379 0 : if(m_nodes[nodeId].childNodeInd[0] == s_invalidIndex)
380 0 : return 0;
381 : return s_numChildren;
382 : }
383 :
384 :
385 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
386 : const size_t* ntree<tree_dim, world_dim, elem_t, common_data_t>::
387 : child_node_ids(size_t nodeId) const
388 : {
389 : assert(nodeId < m_nodes.size());
390 0 : return m_nodes[nodeId].childNodeInd;
391 : }
392 :
393 :
394 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
395 : typename ntree<tree_dim, world_dim, elem_t, common_data_t>::elem_iterator_t
396 : ntree<tree_dim, world_dim, elem_t, common_data_t>::
397 : elems_begin(size_t nodeId) const
398 : {
399 : assert(nodeId < m_nodes.size());
400 0 : if(m_entries.empty())
401 : return elem_iterator_t(NULL, s_invalidIndex);
402 0 : return elem_iterator_t(&m_entries.front(), m_nodes[nodeId].firstEntryInd);
403 : }
404 :
405 :
406 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
407 : typename ntree<tree_dim, world_dim, elem_t, common_data_t>::elem_iterator_t
408 : ntree<tree_dim, world_dim, elem_t, common_data_t>::
409 : elems_end(size_t nodeId) const
410 : {
411 : return elem_iterator_t(NULL, s_invalidIndex);
412 : }
413 :
414 :
415 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
416 : size_t ntree<tree_dim, world_dim, elem_t, common_data_t>::
417 : num_elements(size_t nodeId) const
418 : {
419 : assert(nodeId < m_nodes.size());
420 0 : return m_nodes[nodeId].numEntries;
421 : }
422 :
423 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
424 : size_t ntree<tree_dim, world_dim, elem_t, common_data_t>::
425 : level(size_t nodeId) const
426 : {
427 : assert(nodeId < m_nodes.size());
428 0 : return m_nodes[nodeId].level;
429 : }
430 :
431 : template <int tree_dim, int world_dim, class elem_t, class common_data_t>
432 : const typename ntree<tree_dim, world_dim, elem_t, common_data_t>::box_t&
433 : ntree<tree_dim, world_dim, elem_t, common_data_t>::
434 : bounding_box(size_t nodeId) const
435 : {
436 : assert(nodeId < m_nodes.size());
437 : return m_nodes[nodeId].looseBox;
438 : }
439 :
440 : }// end of namespace
441 :
442 : #endif
|
