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