MOAB: Mesh Oriented datABase  (version 5.2.1)
BSPTree.hpp
Go to the documentation of this file.
00001 /*
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 /**\file BSPTree.hpp
00017  *\author Jason Kraftcheck (kraftche@cae.wisc.edu)
00018  *\date 2008-05-12
00019  */
00020 
00021 #ifndef MOAB_BSP_TREE_HPP
00022 #define MOAB_BSP_TREE_HPP
00023 
00024 #include "moab/Types.hpp"
00025 #include "moab/Interface.hpp"
00026 
00027 #include <math.h>
00028 #include <string>
00029 #include <vector>
00030 
00031 namespace moab
00032 {
00033 
00034 class BSPTreeBoxIter;
00035 class BSPTreeIter;
00036 class Range;
00037 class BSPTreePoly;
00038 
00039 /** \class BSPTree
00040  * \brief BSP tree, for sorting and searching entities spatially
00041  */
00042 class BSPTree
00043 {
00044   private:
00045     Interface* mbInstance;
00046     Tag planeTag, rootTag;
00047     unsigned meshSetFlags;
00048     bool cleanUpTrees;
00049     std::vector< EntityHandle > createdTrees;
00050 
00051     ErrorCode init_tags( const char* tagname = 0 );
00052 
00053   public:
00054     static double epsilon()
00055     {
00056         return 1e-6;
00057     }
00058 
00059     BSPTree( Interface* iface, const char* tagname = 0, unsigned meshset_creation_flags = MESHSET_SET );
00060 
00061     BSPTree( Interface* iface, bool destroy_created_trees, const char* tagname = 0,
00062              unsigned meshset_creation_flags = MESHSET_SET );
00063 
00064     ~BSPTree();
00065 
00066     //! Enumerate split plane directions
00067     enum Axis
00068     {
00069         X = 0,
00070         Y = 1,
00071         Z = 2
00072     };
00073 
00074     /**\brief struct to store a plane
00075      *
00076      * If plane is defined as Ax+By+Cz+D=0, then
00077      * norm={A,B,C} and coeff=-D.
00078      */
00079     struct Plane
00080     {
00081         Plane() : coeff( 0.0 ) {}
00082         Plane( const double n[3], double d ) : coeff( d )
00083         {
00084             norm[0] = n[0];
00085             norm[1] = n[1];
00086             norm[2] = n[2];
00087         }
00088         /** a x + b y + c z + d = 0 */
00089         Plane( double a, double b, double c, double d ) : coeff( d )
00090         {
00091             norm[0] = a;
00092             norm[1] = b;
00093             norm[2] = c;
00094         }
00095         /** Create Y = 1 plane by doing Plane( Y, 1.0 ); */
00096         Plane( Axis normal, double point_on_axis ) : coeff( -point_on_axis )
00097         {
00098             norm[0] = norm[1] = norm[2] = 0;
00099             norm[normal]                = 1.0;
00100         }
00101 
00102         double norm[3];  //!< Unit normal of plane
00103         double coeff;    //!< norm[0]*x + norm[1]*y + norm[2]*z + coeff = 0
00104 
00105         /** Signed distance from point to plane:
00106          *    absolute value is distance from point to plane
00107          *    positive if 'above' plane, negative if 'below'
00108          *  Note: assumes unit-length normal.
00109          */
00110         double signed_distance( const double point[3] ) const
00111         {
00112             return point[0] * norm[0] + point[1] * norm[1] + point[2] * norm[2] + coeff;
00113         }
00114 
00115         /** return true if point is below the plane */
00116         bool below( const double point[3] ) const
00117         {
00118             return signed_distance( point ) <= 0.0;
00119         }
00120 
00121         /** return true if point is above the plane */
00122         bool above( const double point[3] ) const
00123         {
00124             return signed_distance( point ) >= 0.0;
00125         }
00126 
00127         double distance( const double point[3] ) const
00128         {
00129             return fabs( signed_distance( point ) );
00130         }
00131 
00132         /** reverse plane normal */
00133         void flip()
00134         {
00135             norm[0] = -norm[0];
00136             norm[1] = -norm[1];
00137             norm[2] = -norm[2];
00138             coeff   = -coeff;
00139         }
00140 
00141         void set( const double normal[3], const double point[3] )
00142         {
00143             const double dot = normal[0] * point[0] + normal[1] * point[1] + normal[2] * point[2];
00144             *this            = Plane( normal, -dot );
00145         }
00146 
00147         void set( const double pt1[3], const double pt2[3], const double pt3[3] );
00148 
00149         void set( double i, double j, double k, double cff )
00150         {
00151             *this = Plane( i, j, k, cff );
00152         }
00153 
00154         /** Create Y = 1 plane by doing set( Y, 1.0 ); */
00155         void set( Axis normal, double point_on_axis )
00156         {
00157             coeff   = -point_on_axis;
00158             norm[0] = norm[1] = norm[2] = 0;
00159             norm[normal]                = 1.0;
00160         }
00161     };
00162 
00163     //! Get split plane for tree node
00164     ErrorCode get_split_plane( EntityHandle node, Plane& plane )
00165     {
00166         return moab()->tag_get_data( planeTag, &node, 1, &plane );
00167     }
00168 
00169     //! Set split plane for tree node
00170     ErrorCode set_split_plane( EntityHandle node, const Plane& plane );
00171 
00172     //! Get bounding box for entire tree
00173     ErrorCode get_tree_box( EntityHandle root_node, double corner_coords[8][3] );
00174 
00175     //! Get bounding box for entire tree
00176     ErrorCode get_tree_box( EntityHandle root_node, double corner_coords[24] );
00177 
00178     //! Set bounding box for entire tree
00179     ErrorCode set_tree_box( EntityHandle root_node, const double box_min[3], const double box_max[3] );
00180     ErrorCode set_tree_box( EntityHandle root_node, const double corner_coords[8][3] );
00181 
00182     //! Create tree root node
00183     ErrorCode create_tree( const double box_min[3], const double box_max[3], EntityHandle& root_handle );
00184     ErrorCode create_tree( const double corner_coords[8][3], EntityHandle& root_handle );
00185 
00186     //! Create tree root node
00187     ErrorCode create_tree( EntityHandle& root_handle );
00188 
00189     //! Find all tree roots
00190     ErrorCode find_all_trees( Range& results );
00191 
00192     //! Destroy a tree
00193     ErrorCode delete_tree( EntityHandle root_handle );
00194 
00195     Interface* moab()
00196     {
00197         return mbInstance;
00198     }
00199 
00200     //! Get iterator for tree
00201     ErrorCode get_tree_iterator( EntityHandle tree_root, BSPTreeIter& result );
00202 
00203     //! Get iterator at right-most ('last') leaf.
00204     ErrorCode get_tree_end_iterator( EntityHandle tree_root, BSPTreeIter& result );
00205 
00206     //! Split leaf of tree
00207     //! Updates iterator location to point to first new leaf node.
00208     ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane );
00209 
00210     //! Split leaf of tree
00211     //! Updates iterator location to point to first new leaf node.
00212     ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane, EntityHandle& left_child, EntityHandle& right_child );
00213 
00214     //! Split leaf of tree
00215     //! Updates iterator location to point to first new leaf node.
00216     ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane, const Range& left_entities, const Range& right_entities );
00217 
00218     //! Split leaf of tree
00219     //! Updates iterator location to point to first new leaf node.
00220     ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane, const std::vector< EntityHandle >& left_entities,
00221                           const std::vector< EntityHandle >& right_entities );
00222 
00223     //! Merge the leaf pointed to by the current iterator with it's
00224     //! sibling.  If the sibling is not a leaf, multiple merges may
00225     //! be done.
00226     ErrorCode merge_leaf( BSPTreeIter& iter );
00227 
00228     //! Get leaf containing input position.
00229     //!
00230     //! Does not take into account global bounding box of tree.
00231     //! - Therefore there is always one leaf containing the point.
00232     //! - If caller wants to account for global bounding box, then
00233     //!   caller can test against that box and not call this method
00234     //!   at all if the point is outside the box, as there is no leaf
00235     //!   containing the point in that case.
00236     ErrorCode leaf_containing_point( EntityHandle tree_root, const double point[3], EntityHandle& leaf_out );
00237 
00238     //! Get iterator at leaf containing input position.
00239     //!
00240     //! Returns MB_ENTITY_NOT_FOUND if point is not within
00241     //! bounding box of tree.
00242     ErrorCode leaf_containing_point( EntityHandle tree_root, const double xyz[3], BSPTreeIter& result );
00243 };
00244 
00245 /** \class BSPTreeIter
00246  * \brief Iterate over leaves of a BSPTree
00247  */
00248 class BSPTreeIter
00249 {
00250   public:
00251     enum Direction
00252     {
00253         LEFT  = 0,
00254         RIGHT = 1
00255     };
00256 
00257   private:
00258     friend class BSPTree;
00259 
00260     BSPTree* treeTool;
00261 
00262   protected:
00263     std::vector< EntityHandle > mStack;
00264     mutable std::vector< EntityHandle > childVect;
00265 
00266     virtual ErrorCode step_to_first_leaf( Direction direction );
00267 
00268     virtual ErrorCode up();
00269     virtual ErrorCode down( const BSPTree::Plane& plane, Direction direction );
00270 
00271     virtual ErrorCode initialize( BSPTree* tool, EntityHandle root, const double* point = 0 );
00272 
00273   public:
00274     BSPTreeIter() : treeTool( 0 ), childVect( 2 ) {}
00275     virtual ~BSPTreeIter() {}
00276 
00277     BSPTree* tool() const
00278     {
00279         return treeTool;
00280     }
00281 
00282     //! Get handle for current leaf
00283     EntityHandle handle() const
00284     {
00285         return mStack.back();
00286     }
00287 
00288     //! Get depth in tree. root is at depth of 1.
00289     unsigned depth() const
00290     {
00291         return mStack.size();
00292     }
00293 
00294     //! Advance the iterator either left or right in the tree
00295     //! Note:  stepping past the end of the tree will invalidate
00296     //!        the iterator.  It will *not* be work step the
00297     //!        other direction.
00298     virtual ErrorCode step( Direction direction );
00299 
00300     //! Advance to next leaf
00301     //! Returns MB_ENTITY_NOT_FOUND if at end.
00302     //! Note: stepping past the end of the tree will invalidate
00303     //!       the iterator. Calling back() will not work.
00304     ErrorCode step()
00305     {
00306         return step( RIGHT );
00307     }
00308 
00309     //! Move back to previous leaf
00310     //! Returns MB_ENTITY_NOT_FOUND if at beginning.
00311     //! Note: stepping past the start of the tree will invalidate
00312     //!       the iterator. Calling step() will not work.
00313     ErrorCode back()
00314     {
00315         return step( LEFT );
00316     }
00317 
00318     //! Get split plane that separates this node from its immediate sibling.
00319     ErrorCode get_parent_split_plane( BSPTree::Plane& plane ) const;
00320 
00321     //! Get volume of leaf polyhedron
00322     virtual double volume() const;
00323 
00324     //! Find range of overlap between ray and leaf.
00325     //!
00326     //!\param ray_point Coordinates of start point of ray
00327     //!\param ray_vect  Directionion vector for ray such that
00328     //!                 the ray is defined by r(t) = ray_point + t * ray_vect
00329     //!                 for t > 0.
00330     //!\param t_enter   Output: if return value is true, this value
00331     //!                 is the parameter location along the ray at which
00332     //!                 the ray entered the leaf.  If return value is false,
00333     //!                 then this value is undefined.
00334     //!\param t_exit    Output: if return value is true, this value
00335     //!                 is the parameter location along the ray at which
00336     //!                 the ray exited the leaf.  If return value is false,
00337     //!                 then this value is undefined.
00338     //!\return true if ray intersects leaf, false otherwise.
00339     virtual bool intersect_ray( const double ray_point[3], const double ray_vect[3], double& t_enter,
00340                                 double& t_exit ) const;
00341 
00342     //! Return true if thos node and the passed node share the
00343     //! same immediate parent.
00344     bool is_sibling( const BSPTreeIter& other_leaf ) const;
00345 
00346     //! Return true if thos node and the passed node share the
00347     //! same immediate parent.
00348     bool is_sibling( EntityHandle other_leaf ) const;
00349 
00350     //! Returns true if calling step() will advance to the
00351     //! immediate sibling of the current node.  Returns false
00352     //! if current node is root or back() will move to the
00353     //! immediate sibling.
00354     bool sibling_is_forward() const;
00355 
00356     //! Calculate the convex polyhedron bounding this leaf.
00357     virtual ErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const;
00358 };
00359 
00360 /** \class BSPTreeBoxIter
00361  * \brief Iterate over leaves of a BSPTree
00362  */
00363 class BSPTreeBoxIter : public BSPTreeIter
00364 {
00365   private:
00366     double leafCoords[8][3];
00367     struct Corners
00368     {
00369         double coords[4][3];
00370     };
00371     std::vector< Corners > stackData;
00372 
00373   protected:
00374     virtual ErrorCode step_to_first_leaf( Direction direction );
00375 
00376     virtual ErrorCode up();
00377     virtual ErrorCode down( const BSPTree::Plane& plane, Direction direction );
00378 
00379     virtual ErrorCode initialize( BSPTree* tool, EntityHandle root, const double* point = 0 );
00380 
00381   public:
00382     BSPTreeBoxIter() {}
00383     virtual ~BSPTreeBoxIter() {}
00384 
00385     //! Faces of a hex : corner bitmap
00386     enum SideBits
00387     {
00388         B0154 = 0x33,  //!< Face defined by corners {0,1,5,4}: 1st Exodus side
00389         B1265 = 0x66,  //!< Face defined by corners {1,2,6,5}: 2nd Exodus side
00390         B2376 = 0xCC,  //!< Face defined by corners {2,3,7,6}: 3rd Exodus side
00391         B3047 = 0x99,  //!< Face defined by corners {3,0,4,7}: 4th Exodus side
00392         B3210 = 0x0F,  //!< Face defined by corners {3,2,1,0}: 5th Exodus side
00393         B4567 = 0xF0   //!< Face defined by corners {4,5,6,7}: 6th Exodus side
00394     };
00395 
00396     static SideBits side_above_plane( const double hex_coords[8][3], const BSPTree::Plane& plane );
00397 
00398     static SideBits side_on_plane( const double hex_coords[8][3], const BSPTree::Plane& plane );
00399 
00400     static SideBits opposite_face( const SideBits& bits )
00401     {
00402         return ( SideBits )( ( ~bits ) & 0xFF );
00403     }
00404 
00405     static ErrorCode face_corners( const SideBits face, const double hex_corners[8][3], double face_corners_out[4][3] );
00406 
00407     //! Advance the iterator either left or right in the tree
00408     //! Note:  stepping past the end of the tree will invalidate
00409     //!        the iterator.  It will *not* work to subsequently
00410     //!        step the other direction.
00411     virtual ErrorCode step( Direction direction );
00412 
00413     //! Advance to next leaf
00414     //! Returns MB_ENTITY_NOT_FOUND if at end.
00415     //! Note: stepping past the end of the tree will invalidate
00416     //!       the iterator. Calling back() will not work.
00417     ErrorCode step()
00418     {
00419         return BSPTreeIter::step();
00420     }
00421 
00422     //! Move back to previous leaf
00423     //! Returns MB_ENTITY_NOT_FOUND if at beginning.
00424     //! Note: stepping past the start of the tree will invalidate
00425     //!       the iterator. Calling step() will not work.
00426     ErrorCode back()
00427     {
00428         return BSPTreeIter::back();
00429     }
00430 
00431     //! Get coordinates of box corners, in Exodus II hexahedral ordering.
00432     ErrorCode get_box_corners( double coords[8][3] ) const;
00433 
00434     //! Get volume of leaf box
00435     double volume() const;
00436 
00437     //! test if a plane intersects the leaf box
00438     enum XSect
00439     {
00440         MISS   = 0,
00441         SPLIT  = 1,
00442         NONHEX = -1
00443     };
00444     XSect splits( const BSPTree::Plane& plane ) const;
00445 
00446     //! test if a plane intersects the leaf box
00447     bool intersects( const BSPTree::Plane& plane ) const;
00448 
00449     //! Find range of overlap between ray and leaf.
00450     //!
00451     //!\param ray_point Coordinates of start point of ray
00452     //!\param ray_vect  Directionion vector for ray such that
00453     //!                 the ray is defined by r(t) = ray_point + t * ray_vect
00454     //!                 for t > 0.
00455     //!\param t_enter   Output: if return value is true, this value
00456     //!                 is the parameter location along the ray at which
00457     //!                 the ray entered the leaf.  If return value is false,
00458     //!                 then this value is undefined.
00459     //!\param t_exit    Output: if return value is true, this value
00460     //!                 is the parameter location along the ray at which
00461     //!                 the ray exited the leaf.  If return value is false,
00462     //!                 then this value is undefined.
00463     //!\return true if ray intersects leaf, false otherwise.
00464     bool intersect_ray( const double ray_point[3], const double ray_vect[3], double& t_enter, double& t_exit ) const;
00465 
00466     //! Return the side of the box bounding this tree node
00467     //! that is shared with the immediately adjacent sibling
00468     //! (the tree node that shares a common parent node with
00469     //! this node in the binary tree.)
00470     //!
00471     //!\return MB_ENTITY_NOT FOUND if root node.
00472     //!        MB_FAILURE if internal error.
00473     //!        MB_SUCCESS otherwise.
00474     ErrorCode sibling_side( SideBits& side_out ) const;
00475 
00476     //! Get adjacent leaf nodes on indicated side
00477     //!
00478     //!\param side   Face of box for which to retrieve neighbors
00479     //!\param results List to which to append results.  This function does
00480     //!             *not* clear existing values in list.
00481     //!\param epsilon Tolerance on overlap.  A positive value E will
00482     //!              result in nodes that are separated by as much as E
00483     //!              to be considered touching.  A negative value -E will
00484     //!              cause leaves that do not overlap by at least E to be
00485     //!              considered non-overlapping.  Amongst other things,
00486     //!              this value can be used to control whether or not
00487     //!              leaves adjacent at only their edges or corners are
00488     //!              returned.
00489     ErrorCode get_neighbors( SideBits side, std::vector< BSPTreeBoxIter >& results, double epsilon = 0.0 ) const;
00490 
00491     //! Calculate the convex polyhedron bounding this leaf.
00492     ErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const;
00493 };
00494 
00495 }  // namespace moab
00496 
00497 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines