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