Branch data Line data Source code
1 : : /*
2 : : * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 : : * storing and accessing finite element mesh data.
4 : : *
5 : : * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 : : * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 : : * retains certain rights in this software.
8 : : *
9 : : * This library is free software; you can redistribute it and/or
10 : : * modify it under the terms of the GNU Lesser General Public
11 : : * License as published by the Free Software Foundation; either
12 : : * version 2.1 of the License, or (at your option) any later version.
13 : : *
14 : : */
15 : :
16 : : /**\file BSPTree.hpp
17 : : *\author Jason Kraftcheck ([email protected])
18 : : *\date 2008-05-12
19 : : */
20 : :
21 : : #ifndef MOAB_BSP_TREE_HPP
22 : : #define MOAB_BSP_TREE_HPP
23 : :
24 : : #include "moab/Types.hpp"
25 : : #include "moab/Interface.hpp"
26 : :
27 : : #include <math.h>
28 : : #include <string>
29 : : #include <vector>
30 : :
31 : : namespace moab
32 : : {
33 : :
34 : : class BSPTreeBoxIter;
35 : : class BSPTreeIter;
36 : : class Range;
37 : : class BSPTreePoly;
38 : :
39 : : /** \class BSPTree
40 : : * \brief BSP tree, for sorting and searching entities spatially
41 : : */
42 : : class BSPTree
43 : : {
44 : : private:
45 : : Interface* mbInstance;
46 : : Tag planeTag, rootTag;
47 : : unsigned meshSetFlags;
48 : : bool cleanUpTrees;
49 : : std::vector< EntityHandle > createdTrees;
50 : :
51 : : ErrorCode init_tags( const char* tagname = 0 );
52 : :
53 : : public:
54 : 1782 : static double epsilon()
55 : : {
56 : 1782 : return 1e-6;
57 : : }
58 : :
59 : : BSPTree( Interface* iface, const char* tagname = 0, unsigned meshset_creation_flags = MESHSET_SET );
60 : :
61 : : BSPTree( Interface* iface, bool destroy_created_trees, const char* tagname = 0,
62 : : unsigned meshset_creation_flags = MESHSET_SET );
63 : :
64 : : ~BSPTree();
65 : :
66 : : //! Enumerate split plane directions
67 : : enum Axis
68 : : {
69 : : X = 0,
70 : : Y = 1,
71 : : Z = 2
72 : : };
73 : :
74 : : /**\brief struct to store a plane
75 : : *
76 : : * If plane is defined as Ax+By+Cz+D=0, then
77 : : * norm={A,B,C} and coeff=-D.
78 : : */
79 : : struct Plane
80 : : {
81 : 1839 : Plane() : coeff( 0.0 ) {}
82 : 659 : Plane( const double n[3], double d ) : coeff( d )
83 : : {
84 : 659 : norm[0] = n[0];
85 : 659 : norm[1] = n[1];
86 : 659 : norm[2] = n[2];
87 : 659 : }
88 : : /** a x + b y + c z + d = 0 */
89 : 106 : Plane( double a, double b, double c, double d ) : coeff( d )
90 : : {
91 : 106 : norm[0] = a;
92 : 106 : norm[1] = b;
93 : 106 : norm[2] = c;
94 : 106 : }
95 : : /** Create Y = 1 plane by doing Plane( Y, 1.0 ); */
96 : 4 : Plane( Axis normal, double point_on_axis ) : coeff( -point_on_axis )
97 : : {
98 : 4 : norm[0] = norm[1] = norm[2] = 0;
99 : 4 : norm[normal] = 1.0;
100 : 4 : }
101 : :
102 : : double norm[3]; //!< Unit normal of plane
103 : : double coeff; //!< norm[0]*x + norm[1]*y + norm[2]*z + coeff = 0
104 : :
105 : : /** Signed distance from point to plane:
106 : : * absolute value is distance from point to plane
107 : : * positive if 'above' plane, negative if 'below'
108 : : * Note: assumes unit-length normal.
109 : : */
110 : 5360 : double signed_distance( const double point[3] ) const
111 : : {
112 : 5360 : return point[0] * norm[0] + point[1] * norm[1] + point[2] * norm[2] + coeff;
113 : : }
114 : :
115 : : /** return true if point is below the plane */
116 : : bool below( const double point[3] ) const
117 : : {
118 : : return signed_distance( point ) <= 0.0;
119 : : }
120 : :
121 : : /** return true if point is above the plane */
122 : 3439 : bool above( const double point[3] ) const
123 : : {
124 : 3439 : return signed_distance( point ) >= 0.0;
125 : : }
126 : :
127 : 1683 : double distance( const double point[3] ) const
128 : : {
129 : 1683 : return fabs( signed_distance( point ) );
130 : : }
131 : :
132 : : /** reverse plane normal */
133 : 180 : void flip()
134 : : {
135 : 180 : norm[0] = -norm[0];
136 : 180 : norm[1] = -norm[1];
137 : 180 : norm[2] = -norm[2];
138 : 180 : coeff = -coeff;
139 : 180 : }
140 : :
141 : 4 : void set( const double normal[3], const double point[3] )
142 : : {
143 : 4 : const double dot = normal[0] * point[0] + normal[1] * point[1] + normal[2] * point[2];
144 : 4 : *this = Plane( normal, -dot );
145 : 4 : }
146 : :
147 : : void set( const double pt1[3], const double pt2[3], const double pt3[3] );
148 : :
149 : 5 : void set( double i, double j, double k, double cff )
150 : : {
151 : 5 : *this = Plane( i, j, k, cff );
152 : 5 : }
153 : :
154 : : /** Create Y = 1 plane by doing set( Y, 1.0 ); */
155 : : void set( Axis normal, double point_on_axis )
156 : : {
157 : : coeff = -point_on_axis;
158 : : norm[0] = norm[1] = norm[2] = 0;
159 : : norm[normal] = 1.0;
160 : : }
161 : : };
162 : :
163 : : //! Get split plane for tree node
164 : 800 : ErrorCode get_split_plane( EntityHandle node, Plane& plane )
165 : : {
166 : 800 : return moab()->tag_get_data( planeTag, &node, 1, &plane );
167 : : }
168 : :
169 : : //! Set split plane for tree node
170 : : ErrorCode set_split_plane( EntityHandle node, const Plane& plane );
171 : :
172 : : //! Get bounding box for entire tree
173 : : ErrorCode get_tree_box( EntityHandle root_node, double corner_coords[8][3] );
174 : :
175 : : //! Get bounding box for entire tree
176 : : ErrorCode get_tree_box( EntityHandle root_node, double corner_coords[24] );
177 : :
178 : : //! Set bounding box for entire tree
179 : : ErrorCode set_tree_box( EntityHandle root_node, const double box_min[3], const double box_max[3] );
180 : : ErrorCode set_tree_box( EntityHandle root_node, const double corner_coords[8][3] );
181 : :
182 : : //! Create tree root node
183 : : ErrorCode create_tree( const double box_min[3], const double box_max[3], EntityHandle& root_handle );
184 : : ErrorCode create_tree( const double corner_coords[8][3], EntityHandle& root_handle );
185 : :
186 : : //! Create tree root node
187 : : ErrorCode create_tree( EntityHandle& root_handle );
188 : :
189 : : //! Find all tree roots
190 : : ErrorCode find_all_trees( Range& results );
191 : :
192 : : //! Destroy a tree
193 : : ErrorCode delete_tree( EntityHandle root_handle );
194 : :
195 : 2921 : Interface* moab()
196 : : {
197 : 2921 : return mbInstance;
198 : : }
199 : :
200 : : //! Get iterator for tree
201 : : ErrorCode get_tree_iterator( EntityHandle tree_root, BSPTreeIter& result );
202 : :
203 : : //! Get iterator at right-most ('last') leaf.
204 : : ErrorCode get_tree_end_iterator( EntityHandle tree_root, BSPTreeIter& result );
205 : :
206 : : //! Split leaf of tree
207 : : //! Updates iterator location to point to first new leaf node.
208 : : ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane );
209 : :
210 : : //! Split leaf of tree
211 : : //! Updates iterator location to point to first new leaf node.
212 : : ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane, EntityHandle& left_child, EntityHandle& right_child );
213 : :
214 : : //! Split leaf of tree
215 : : //! Updates iterator location to point to first new leaf node.
216 : : ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane, const Range& left_entities, const Range& right_entities );
217 : :
218 : : //! Split leaf of tree
219 : : //! Updates iterator location to point to first new leaf node.
220 : : ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane, const std::vector< EntityHandle >& left_entities,
221 : : const std::vector< EntityHandle >& right_entities );
222 : :
223 : : //! Merge the leaf pointed to by the current iterator with it's
224 : : //! sibling. If the sibling is not a leaf, multiple merges may
225 : : //! be done.
226 : : ErrorCode merge_leaf( BSPTreeIter& iter );
227 : :
228 : : //! Get leaf containing input position.
229 : : //!
230 : : //! Does not take into account global bounding box of tree.
231 : : //! - Therefore there is always one leaf containing the point.
232 : : //! - If caller wants to account for global bounding box, then
233 : : //! caller can test against that box and not call this method
234 : : //! at all if the point is outside the box, as there is no leaf
235 : : //! containing the point in that case.
236 : : ErrorCode leaf_containing_point( EntityHandle tree_root, const double point[3], EntityHandle& leaf_out );
237 : :
238 : : //! Get iterator at leaf containing input position.
239 : : //!
240 : : //! Returns MB_ENTITY_NOT_FOUND if point is not within
241 : : //! bounding box of tree.
242 : : ErrorCode leaf_containing_point( EntityHandle tree_root, const double xyz[3], BSPTreeIter& result );
243 : : };
244 : :
245 : : /** \class BSPTreeIter
246 : : * \brief Iterate over leaves of a BSPTree
247 : : */
248 [ + - ]: 168 : class BSPTreeIter
249 : : {
250 : : public:
251 : : enum Direction
252 : : {
253 : : LEFT = 0,
254 : : RIGHT = 1
255 : : };
256 : :
257 : : private:
258 : : friend class BSPTree;
259 : :
260 : : BSPTree* treeTool;
261 : :
262 : : protected:
263 : : std::vector< EntityHandle > mStack;
264 : : mutable std::vector< EntityHandle > childVect;
265 : :
266 : : virtual ErrorCode step_to_first_leaf( Direction direction );
267 : :
268 : : virtual ErrorCode up();
269 : : virtual ErrorCode down( const BSPTree::Plane& plane, Direction direction );
270 : :
271 : : virtual ErrorCode initialize( BSPTree* tool, EntityHandle root, const double* point = 0 );
272 : :
273 : : public:
274 [ + - ]: 27 : BSPTreeIter() : treeTool( 0 ), childVect( 2 ) {}
275 [ - + ]: 342 : virtual ~BSPTreeIter() {}
276 : :
277 : 1521 : BSPTree* tool() const
278 : : {
279 : 1521 : return treeTool;
280 : : }
281 : :
282 : : //! Get handle for current leaf
283 : 1190 : EntityHandle handle() const
284 : : {
285 : 1190 : return mStack.back();
286 : : }
287 : :
288 : : //! Get depth in tree. root is at depth of 1.
289 : 46 : unsigned depth() const
290 : : {
291 : 46 : return mStack.size();
292 : : }
293 : :
294 : : //! Advance the iterator either left or right in the tree
295 : : //! Note: stepping past the end of the tree will invalidate
296 : : //! the iterator. It will *not* be work step the
297 : : //! other direction.
298 : : virtual ErrorCode step( Direction direction );
299 : :
300 : : //! Advance to next leaf
301 : : //! Returns MB_ENTITY_NOT_FOUND if at end.
302 : : //! Note: stepping past the end of the tree will invalidate
303 : : //! the iterator. Calling back() will not work.
304 : 104 : ErrorCode step()
305 : : {
306 : 104 : return step( RIGHT );
307 : : }
308 : :
309 : : //! Move back to previous leaf
310 : : //! Returns MB_ENTITY_NOT_FOUND if at beginning.
311 : : //! Note: stepping past the start of the tree will invalidate
312 : : //! the iterator. Calling step() will not work.
313 : 6 : ErrorCode back()
314 : : {
315 : 6 : return step( LEFT );
316 : : }
317 : :
318 : : //! Get split plane that separates this node from its immediate sibling.
319 : : ErrorCode get_parent_split_plane( BSPTree::Plane& plane ) const;
320 : :
321 : : //! Get volume of leaf polyhedron
322 : : virtual double volume() const;
323 : :
324 : : //! Find range of overlap between ray and leaf.
325 : : //!
326 : : //!\param ray_point Coordinates of start point of ray
327 : : //!\param ray_vect Directionion vector for ray such that
328 : : //! the ray is defined by r(t) = ray_point + t * ray_vect
329 : : //! for t > 0.
330 : : //!\param t_enter Output: if return value is true, this value
331 : : //! is the parameter location along the ray at which
332 : : //! the ray entered the leaf. If return value is false,
333 : : //! then this value is undefined.
334 : : //!\param t_exit Output: if return value is true, this value
335 : : //! is the parameter location along the ray at which
336 : : //! the ray exited the leaf. If return value is false,
337 : : //! then this value is undefined.
338 : : //!\return true if ray intersects leaf, false otherwise.
339 : : virtual bool intersect_ray( const double ray_point[3], const double ray_vect[3], double& t_enter,
340 : : double& t_exit ) const;
341 : :
342 : : //! Return true if thos node and the passed node share the
343 : : //! same immediate parent.
344 : : bool is_sibling( const BSPTreeIter& other_leaf ) const;
345 : :
346 : : //! Return true if thos node and the passed node share the
347 : : //! same immediate parent.
348 : : bool is_sibling( EntityHandle other_leaf ) const;
349 : :
350 : : //! Returns true if calling step() will advance to the
351 : : //! immediate sibling of the current node. Returns false
352 : : //! if current node is root or back() will move to the
353 : : //! immediate sibling.
354 : : bool sibling_is_forward() const;
355 : :
356 : : //! Calculate the convex polyhedron bounding this leaf.
357 : : virtual ErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const;
358 : : };
359 : :
360 : : /** \class BSPTreeBoxIter
361 : : * \brief Iterate over leaves of a BSPTree
362 : : */
363 [ + - ]: 168 : class BSPTreeBoxIter : public BSPTreeIter
364 : : {
365 : : private:
366 : : double leafCoords[8][3];
367 : : struct Corners
368 : : {
369 : : double coords[4][3];
370 : : };
371 : : std::vector< Corners > stackData;
372 : :
373 : : protected:
374 : : virtual ErrorCode step_to_first_leaf( Direction direction );
375 : :
376 : : virtual ErrorCode up();
377 : : virtual ErrorCode down( const BSPTree::Plane& plane, Direction direction );
378 : :
379 : : virtual ErrorCode initialize( BSPTree* tool, EntityHandle root, const double* point = 0 );
380 : :
381 : : public:
382 [ + - ]: 10 : BSPTreeBoxIter() {}
383 [ - + ]: 308 : virtual ~BSPTreeBoxIter() {}
384 : :
385 : : //! Faces of a hex : corner bitmap
386 : : enum SideBits
387 : : {
388 : : B0154 = 0x33, //!< Face defined by corners {0,1,5,4}: 1st Exodus side
389 : : B1265 = 0x66, //!< Face defined by corners {1,2,6,5}: 2nd Exodus side
390 : : B2376 = 0xCC, //!< Face defined by corners {2,3,7,6}: 3rd Exodus side
391 : : B3047 = 0x99, //!< Face defined by corners {3,0,4,7}: 4th Exodus side
392 : : B3210 = 0x0F, //!< Face defined by corners {3,2,1,0}: 5th Exodus side
393 : : B4567 = 0xF0 //!< Face defined by corners {4,5,6,7}: 6th Exodus side
394 : : };
395 : :
396 : : static SideBits side_above_plane( const double hex_coords[8][3], const BSPTree::Plane& plane );
397 : :
398 : : static SideBits side_on_plane( const double hex_coords[8][3], const BSPTree::Plane& plane );
399 : :
400 : 77 : static SideBits opposite_face( const SideBits& bits )
401 : : {
402 : 77 : return ( SideBits )( ( ~bits ) & 0xFF );
403 : : }
404 : :
405 : : static ErrorCode face_corners( const SideBits face, const double hex_corners[8][3], double face_corners_out[4][3] );
406 : :
407 : : //! Advance the iterator either left or right in the tree
408 : : //! Note: stepping past the end of the tree will invalidate
409 : : //! the iterator. It will *not* work to subsequently
410 : : //! step the other direction.
411 : : virtual ErrorCode step( Direction direction );
412 : :
413 : : //! Advance to next leaf
414 : : //! Returns MB_ENTITY_NOT_FOUND if at end.
415 : : //! Note: stepping past the end of the tree will invalidate
416 : : //! the iterator. Calling back() will not work.
417 : 30 : ErrorCode step()
418 : : {
419 : 30 : return BSPTreeIter::step();
420 : : }
421 : :
422 : : //! Move back to previous leaf
423 : : //! Returns MB_ENTITY_NOT_FOUND if at beginning.
424 : : //! Note: stepping past the start of the tree will invalidate
425 : : //! the iterator. Calling step() will not work.
426 : 2 : ErrorCode back()
427 : : {
428 : 2 : return BSPTreeIter::back();
429 : : }
430 : :
431 : : //! Get coordinates of box corners, in Exodus II hexahedral ordering.
432 : : ErrorCode get_box_corners( double coords[8][3] ) const;
433 : :
434 : : //! Get volume of leaf box
435 : : double volume() const;
436 : :
437 : : //! test if a plane intersects the leaf box
438 : : enum XSect
439 : : {
440 : : MISS = 0,
441 : : SPLIT = 1,
442 : : NONHEX = -1
443 : : };
444 : : XSect splits( const BSPTree::Plane& plane ) const;
445 : :
446 : : //! test if a plane intersects the leaf box
447 : : bool intersects( const BSPTree::Plane& plane ) const;
448 : :
449 : : //! Find range of overlap between ray and leaf.
450 : : //!
451 : : //!\param ray_point Coordinates of start point of ray
452 : : //!\param ray_vect Directionion vector for ray such that
453 : : //! the ray is defined by r(t) = ray_point + t * ray_vect
454 : : //! for t > 0.
455 : : //!\param t_enter Output: if return value is true, this value
456 : : //! is the parameter location along the ray at which
457 : : //! the ray entered the leaf. If return value is false,
458 : : //! then this value is undefined.
459 : : //!\param t_exit Output: if return value is true, this value
460 : : //! is the parameter location along the ray at which
461 : : //! the ray exited the leaf. If return value is false,
462 : : //! then this value is undefined.
463 : : //!\return true if ray intersects leaf, false otherwise.
464 : : bool intersect_ray( const double ray_point[3], const double ray_vect[3], double& t_enter, double& t_exit ) const;
465 : :
466 : : //! Return the side of the box bounding this tree node
467 : : //! that is shared with the immediately adjacent sibling
468 : : //! (the tree node that shares a common parent node with
469 : : //! this node in the binary tree.)
470 : : //!
471 : : //!\return MB_ENTITY_NOT FOUND if root node.
472 : : //! MB_FAILURE if internal error.
473 : : //! MB_SUCCESS otherwise.
474 : : ErrorCode sibling_side( SideBits& side_out ) const;
475 : :
476 : : //! Get adjacent leaf nodes on indicated side
477 : : //!
478 : : //!\param side Face of box for which to retrieve neighbors
479 : : //!\param results List to which to append results. This function does
480 : : //! *not* clear existing values in list.
481 : : //!\param epsilon Tolerance on overlap. A positive value E will
482 : : //! result in nodes that are separated by as much as E
483 : : //! to be considered touching. A negative value -E will
484 : : //! cause leaves that do not overlap by at least E to be
485 : : //! considered non-overlapping. Amongst other things,
486 : : //! this value can be used to control whether or not
487 : : //! leaves adjacent at only their edges or corners are
488 : : //! returned.
489 : : ErrorCode get_neighbors( SideBits side, std::vector< BSPTreeBoxIter >& results, double epsilon = 0.0 ) const;
490 : :
491 : : //! Calculate the convex polyhedron bounding this leaf.
492 : : ErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const;
493 : : };
494 : :
495 : : } // namespace moab
496 : :
497 : : #endif
|
