Line data Source code
1 : /*
2 : * Copyright (c) 2010-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 : #include <cassert>
34 : #include <vector>
35 : #include "octree.h"
36 :
37 : using namespace std;
38 :
39 : namespace ug{
40 : namespace node_tree
41 : {
42 :
43 : ////////////////////////////////////////////////////////////////////////
44 : /** calculates the bounding box around a set of points*/
45 : static void
46 0 : CalculateBoundingBox(vector3& boxMinOut, vector3& boxMaxOut,
47 : const vector3* points, size_t numPoints)
48 : {
49 0 : if(numPoints < 1){
50 : boxMinOut = boxMaxOut = vector3(0, 0, 0);
51 0 : return;
52 : }
53 :
54 : boxMinOut = boxMaxOut = points[0];
55 :
56 0 : for(size_t i = 0; i < numPoints; ++i){
57 0 : for(size_t j = 0; j < 3; ++j)
58 : {
59 0 : const number& coord = points[i][j];
60 0 : if(coord < boxMinOut[j])
61 0 : boxMinOut[j] = coord;
62 0 : else if(coord > boxMaxOut[j])
63 0 : boxMaxOut[j] = coord;
64 : }
65 : }
66 : }
67 :
68 : ////////////////////////////////////////////////////////////////////////
69 : /** calculates the bounding box around a set of elements which are
70 : * given by a set of indices referencing a set of points. This method
71 : * is particularily useful if the elements only reference a subset
72 : * of the given points. If not the bounding box should be calculated
73 : * directly on the set of points.
74 : */
75 : static void
76 0 : CalculateBoundingBox(vector3& boxMinOut, vector3& boxMaxOut,
77 : const vector3* points,
78 : const int* inds, size_t numInds)
79 : {
80 0 : if(numInds < 1){
81 : boxMinOut = boxMaxOut = vector3(0, 0, 0);
82 0 : return;
83 : }
84 :
85 0 : boxMinOut = boxMaxOut = points[inds[0]];
86 :
87 0 : for(size_t i = 0; i < numInds; ++i){
88 0 : for(size_t j = 0; j < 3; ++j)
89 : {
90 0 : const number& coord = points[inds[i]][j];
91 0 : if(coord < boxMinOut[j])
92 0 : boxMinOut[j] = coord;
93 0 : else if(coord > boxMaxOut[j])
94 0 : boxMaxOut[j] = coord;
95 : }
96 : }
97 : }
98 :
99 : ////////////////////////////////////////////////////////////////////////
100 : /** Increases the size of the bounding box by scaling around its center*/
101 0 : static void GrowBox(vector3& minInOut, vector3& maxInOut, number scaleFac)
102 : {
103 : vector3 halfDiag;
104 : VecSubtract(halfDiag, maxInOut, minInOut);
105 : VecScale(halfDiag, halfDiag, 0.5);
106 : VecScale(halfDiag, halfDiag, scaleFac);
107 :
108 : vector3 center;
109 : VecAdd(center, minInOut, maxInOut);
110 : VecScale(center, center, 0.5);
111 :
112 : VecAdd(maxInOut, center, halfDiag);
113 : VecSubtract(minInOut, center, halfDiag);
114 0 : }
115 :
116 : ////////////////////////////////////////////////////////////////////////
117 : // CreateOctree
118 : SPCollisionTreeRootNode
119 0 : CreateOctree(vector3* points, size_t numPoints,
120 : int* elemInds, size_t numElemInds, int numIndsPerElem,
121 : CollisionElementID* elemIDs, int maxDepth,
122 : int elemThreshold, bool bLoose)
123 : {
124 : // if there are no elements, we'll return an empty node
125 0 : if(numElemInds == 0)
126 : return SPCollisionTreeRootNode(NULL);
127 :
128 : // if the number of indices per elements is not right, we'll return, too.
129 0 : if(numIndsPerElem < 2 || numIndsPerElem > 3)
130 : return SPCollisionTreeRootNode(NULL);
131 :
132 : // calculate the bounding-box of the tree
133 : vector3 boxMin, boxMax;
134 0 : CalculateBoundingBox(boxMin, boxMax, points, numPoints);
135 0 : GrowBox(boxMin, boxMax, 1.001);
136 :
137 : // create the root-node
138 0 : SPCollisionTreeRootNode spRootNode = CollisionTreeRootNode::create();
139 :
140 : // add the given points to the root-node
141 0 : spRootNode->add_points(points, numPoints);
142 :
143 : // set the bounding-box of the root-node
144 0 : spRootNode->set_box(boxMin, boxMax);
145 :
146 : // create the sub-trees
147 0 : CreateSubOctrees(spRootNode.get(), points, numPoints, elemInds,
148 : numElemInds, numIndsPerElem, elemIDs, maxDepth,
149 : elemThreshold, bLoose);
150 :
151 : // done. return the root-node
152 : return spRootNode;
153 : }
154 :
155 : ////////////////////////////////////////////////////////////////////////
156 : // CreateSubOctrees
157 0 : void CreateSubOctrees(BoxedGroupNode* parentNode,
158 : vector3* points, size_t numPoints,
159 : int* elemInds, size_t numElemInds, int numIndsPerElem,
160 : CollisionElementID* elemIDs, int maxDepth,
161 : int elemThreshold, bool bLoose)
162 : {
163 : // the maximal number of indices per element.
164 : // has to be adjusted if new element-types are supported.
165 : assert(numIndsPerElem > 1 && numIndsPerElem <= 3 &&
166 : "unsupported number of indices per element.");
167 :
168 : // the number of elements
169 0 : int numElems = numElemInds / numIndsPerElem;
170 :
171 : // check if subnodes have to be created
172 0 : if((maxDepth > 0) && (numElems > elemThreshold))
173 : {
174 : // yes, we have to create subnodes.
175 : vector3 boxMin, boxMax;
176 0 : boxMin = parentNode->min_corner();
177 0 : boxMax = parentNode->max_corner();
178 :
179 : vector3 boxSize;
180 0 : boxSize.x() = boxMax.x() - boxMin.x();
181 0 : boxSize.y() = boxMax.y() - boxMin.y();
182 0 : boxSize.z() = boxMax.z() - boxMin.z();
183 :
184 : // we'll store index-lists of subtrees in those vectors
185 : vector<int> vNewElems;
186 : vector<CollisionElementID> vNewIDs;
187 : //TODO: check whether vNewElems.reserve amd vNewIDs.reserve would result in a speedup.
188 :
189 : // if a loose tree shall be generated, we'll have to record
190 : // whether an edge has already been assigned or not.
191 : vector<char> elemAssigned;
192 0 : if(bLoose)
193 0 : elemAssigned.resize(numElems, 0);
194 :
195 : // iterate through each sub-tree
196 0 : for(size_t subNodeInd = 0; subNodeInd < 8; ++subNodeInd)
197 : {
198 : // clear arrays
199 : vNewElems.clear();
200 : vNewIDs.clear();
201 :
202 : // create the box of the subnode
203 : vector3 subBoxMin, subBoxMax;
204 :
205 : // compute bounding box from ParentBox and index
206 0 : subBoxMin.x() = boxMin.x() + .5 * boxSize.x() * (subNodeInd % 2);
207 0 : subBoxMin.y() = boxMin.y() + .5 * boxSize.y() * (int)(subNodeInd / 4);
208 0 : subBoxMin.z() = boxMin.z() + .5 * boxSize.z() * (int)((subNodeInd % 4) / 2);
209 0 : subBoxMax.x() = subBoxMin.x() + .5 * boxSize.x();
210 0 : subBoxMax.y() = subBoxMin.y() + .5 * boxSize.y();
211 0 : subBoxMax.z() = subBoxMin.z() + .5 * boxSize.z();
212 :
213 : // find the elements that are inside the box.
214 : // different approaches have to be taken, depending on
215 : // whether a loose tree shall be created or not.
216 0 : if(bLoose){
217 0 : GrowBox(subBoxMin, subBoxMax, 1.001);
218 : // create a loose tree.
219 : // check for each element whether the center lies in the box or not.
220 0 : for(size_t i = 0; i < numElemInds; i += numIndsPerElem)
221 : {
222 0 : size_t elemIndex = i / numIndsPerElem;
223 :
224 0 : if(!elemAssigned[elemIndex])
225 : {
226 : // calculate the center
227 : vector3 center(0, 0, 0);
228 0 : for(int j = 0; j < numIndsPerElem; ++j){
229 0 : VecAdd(center, center, points[elemInds[i + j]]);
230 : }
231 :
232 0 : VecScale(center, center, 1. / numIndsPerElem);
233 :
234 : // check whether the center lies in the bounding box
235 0 : if(BoxBoundProbe(center, subBoxMin, subBoxMax))
236 : {
237 : // it does. Mark the element as assigned and
238 : // add it to the new-lists.
239 0 : elemAssigned[elemIndex] = 1;
240 :
241 0 : for(int j = 0; j < numIndsPerElem; ++j)
242 0 : vNewElems.push_back(elemInds[i + j]);
243 :
244 0 : if(elemIDs)
245 0 : vNewIDs.push_back(elemIDs[elemIndex]);
246 : }
247 : }
248 : }
249 :
250 : // all elements for this sub-tree have been collected.
251 : // Since it is a loose tree, we'll have to readjust the box
252 0 : CalculateBoundingBox(subBoxMin, subBoxMax, points,
253 : &vNewElems.front(), vNewElems.size());
254 :
255 0 : GrowBox(subBoxMin, subBoxMax, 1.001);
256 : }
257 : else{
258 : //TODO: check whether moving the element-iteration into the switch would
259 : // result in a noticeable speedup (I don't think so!).
260 :
261 : // the tree is not loose. We have to find all elements that
262 : // intersect the box.
263 0 : for(size_t i = 0; i < numElemInds; i += numIndsPerElem)
264 : {
265 : // check whether elements are lines or triangles
266 : bool intersects = false;
267 0 : switch(numIndsPerElem)
268 : {
269 0 : case 2: // elements are line-segments
270 0 : intersects = LineBoxIntersection(points[elemInds[i]],
271 0 : points[elemInds[i + 1]],
272 : subBoxMin, subBoxMax);
273 0 : break;
274 :
275 0 : case 3: // elements are triangles
276 0 : intersects = TriangleBoxIntersection(
277 0 : points[elemInds[i]], points[elemInds[i + 1]],
278 0 : points[elemInds[i + 2]], subBoxMin, subBoxMax);
279 : break;
280 : default:
281 : intersects = false;
282 : break;
283 : }
284 :
285 : // if we found an intersection, the element has to be assigned
286 : // to the subtree, which is associated with the box.
287 0 : if(intersects){
288 0 : for(int j = 0; j < numIndsPerElem; ++j)
289 0 : vNewElems.push_back(elemInds[i + j]);
290 :
291 0 : if(elemIDs)
292 0 : vNewIDs.push_back(elemIDs[i / numIndsPerElem]);
293 : }
294 : }
295 : }
296 :
297 : // all elements that go into this subtree have been found.
298 : // now create the subtree
299 0 : SPBoxedGroupNode spSubNode = BoxedGroupNode::create();
300 0 : spSubNode->set_box(subBoxMin, subBoxMax);
301 0 : parentNode->add_child(spSubNode);
302 :
303 0 : if(elemIDs)
304 0 : CreateSubOctrees(spSubNode.get(), points, numPoints,
305 : &vNewElems.front(), vNewElems.size(), numIndsPerElem,
306 : &vNewIDs.front(), maxDepth - 1, elemThreshold, bLoose);
307 : else
308 0 : CreateSubOctrees(spSubNode.get(), points, numPoints,
309 : &vNewElems.front(), vNewElems.size(), numIndsPerElem,
310 : NULL, maxDepth - 1, elemThreshold, bLoose);
311 : }
312 0 : }
313 : else{
314 : // there are less elements thatn elemThreshold or maxDepth is reached.
315 : // We will now create the tree-node that holds the elements.
316 : // If there are no elements, we can return immediatly.
317 0 : if(numElemInds == 0)
318 : return;
319 :
320 : // We have to distinguish between the different element-types
321 0 : switch(numIndsPerElem){
322 0 : case 2:// edges
323 : {
324 0 : SPCollisionEdgesNode spNode = CollisionEdgesNode::create();
325 0 : if(elemIDs)
326 0 : spNode->add_edges(elemInds, elemIDs, numElemInds / 2);
327 : else
328 0 : spNode->add_edges(elemInds, numElemInds / 2);
329 :
330 0 : parentNode->add_child(spNode);
331 : }
332 0 : break;
333 :
334 0 : case 3:// triangles
335 : {
336 0 : SPCollisionTrianglesNode spNode = CollisionTrianglesNode::create();
337 0 : if(elemIDs)
338 0 : spNode->add_triangles(elemInds, elemIDs, numElemInds / 3);
339 : else
340 0 : spNode->add_triangles(elemInds, numElemInds / 3);
341 :
342 0 : parentNode->add_child(spNode);
343 : }
344 0 : break;
345 :
346 : default:
347 : break;
348 : }
349 : }
350 : }
351 :
352 :
353 : }// end of namesapce node_tree
354 : }// end of namespace ug
355 :
|
