moab
|
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 Coroporation, 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 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 class BSPTreeBoxIter; 00034 class BSPTreeIter; 00035 class Range; 00036 class BSPTreePoly; 00037 00041 class BSPTree 00042 { 00043 private: 00044 Interface* mbInstance; 00045 Tag planeTag, rootTag; 00046 unsigned meshSetFlags; 00047 bool cleanUpTrees; 00048 std::vector<EntityHandle> createdTrees; 00049 00050 ErrorCode init_tags( const char* tagname = 0 ); 00051 00052 public: 00053 00054 static double epsilon() { return 1e-6; } 00055 00056 BSPTree( Interface* iface, 00057 const char* tagname = 0, 00058 unsigned meshset_creation_flags = MESHSET_SET ); 00059 00060 BSPTree( Interface* iface, 00061 bool destroy_created_trees, 00062 const char* tagname = 0, 00063 unsigned meshset_creation_flags = MESHSET_SET ); 00064 00065 ~BSPTree(); 00066 00068 enum Axis { X = 0, Y = 1, Z = 2 }; 00069 00075 struct Plane { 00076 Plane() {} 00077 Plane( const double n[3], double d ) : coeff(d) 00078 { norm[0] = n[0]; norm[1] = n[1]; norm[2] = n[2]; } 00080 Plane( double a, double b, double c, double d ) : coeff(d) 00081 { norm[0] = a; norm[1] = b; norm[2] = c; } 00083 Plane( Axis normal, double point_on_axis ) : coeff(-point_on_axis) 00084 { norm[0] = norm[1] = norm[2] = 0; norm[normal] = 1.0; } 00085 00086 double norm[3]; 00087 double coeff; 00088 00094 double signed_distance( const double point[3] ) const 00095 { return point[0]*norm[0]+point[1]*norm[1]+point[2]*norm[2] + coeff; } 00096 00098 bool below( const double point[3] ) const 00099 { return signed_distance(point) <= 0.0; } 00100 00102 bool above( const double point[3] ) const 00103 { return signed_distance(point) >= 0.0; } 00104 00105 double distance( const double point[3] ) const 00106 { return fabs( signed_distance(point) ); } 00107 00109 void flip() 00110 { norm[0] = -norm[0]; norm[1] = -norm[1]; norm[2] = -norm[2]; coeff = -coeff; } 00111 00112 void set( const double normal[3], const double point[3] ) 00113 { 00114 const double dot = normal[0]*point[0] + normal[1]*point[1] + normal[2]*point[2]; 00115 *this = Plane( normal, -dot ); 00116 } 00117 00118 void set( const double pt1[3], const double pt2[3], const double pt3[3] ); 00119 00120 void set( double i, double j, double k, double cff ) 00121 { *this = Plane( i, j, k, cff ); } 00122 00124 void set( Axis normal, double point_on_axis ) 00125 { 00126 coeff = -point_on_axis; 00127 norm[0] = norm[1] = norm[2] = 0; 00128 norm[normal] = 1.0; 00129 } 00130 }; 00131 00133 ErrorCode get_split_plane( EntityHandle node, Plane& plane ) 00134 { return moab()->tag_get_data( planeTag, &node, 1, &plane ); } 00135 00137 ErrorCode set_split_plane( EntityHandle node, const Plane& plane ); 00138 00140 ErrorCode get_tree_box( EntityHandle root_node, 00141 double corner_coords[8][3] ); 00142 00144 ErrorCode get_tree_box( EntityHandle root_node, 00145 double corner_coords[24] ); 00146 00148 ErrorCode set_tree_box( EntityHandle root_node, 00149 const double box_min[3], 00150 const double box_max[3] ); 00151 ErrorCode set_tree_box( EntityHandle root_node, 00152 const double corner_coords[8][3] ); 00153 00155 ErrorCode create_tree( const double box_min[3], 00156 const double box_max[3], 00157 EntityHandle& root_handle ); 00158 ErrorCode create_tree( const double corner_coords[8][3], 00159 EntityHandle& root_handle ); 00160 00162 ErrorCode create_tree( EntityHandle& root_handle ); 00163 00165 ErrorCode find_all_trees( Range& results ); 00166 00168 ErrorCode delete_tree( EntityHandle root_handle ); 00169 00170 Interface* moab() { return mbInstance; } 00171 00173 ErrorCode get_tree_iterator( EntityHandle tree_root, 00174 BSPTreeIter& result ); 00175 00177 ErrorCode get_tree_end_iterator( EntityHandle tree_root, 00178 BSPTreeIter& result ); 00179 00182 ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane ); 00183 00186 ErrorCode split_leaf( BSPTreeIter& leaf, 00187 Plane plane, 00188 EntityHandle& left_child, 00189 EntityHandle& right_child ); 00190 00193 ErrorCode split_leaf( BSPTreeIter& leaf, 00194 Plane plane, 00195 const Range& left_entities, 00196 const Range& right_entities ); 00197 00200 ErrorCode split_leaf( BSPTreeIter& leaf, 00201 Plane plane, 00202 const std::vector<EntityHandle>& left_entities, 00203 const std::vector<EntityHandle>& right_entities ); 00204 00208 ErrorCode merge_leaf( BSPTreeIter& iter ); 00209 00218 ErrorCode leaf_containing_point( EntityHandle tree_root, 00219 const double point[3], 00220 EntityHandle& leaf_out ); 00221 00226 ErrorCode leaf_containing_point( EntityHandle tree_root, 00227 const double xyz[3], 00228 BSPTreeIter& result ); 00229 }; 00230 00234 class BSPTreeIter 00235 { 00236 public: 00237 00238 enum Direction { LEFT = 0, RIGHT = 1 }; 00239 00240 private: 00241 friend class BSPTree; 00242 00243 BSPTree* treeTool; 00244 00245 protected: 00246 std::vector<EntityHandle> mStack; 00247 mutable std::vector<EntityHandle> childVect; 00248 00249 virtual ErrorCode step_to_first_leaf( Direction direction ); 00250 00251 virtual ErrorCode up(); 00252 virtual ErrorCode down( const BSPTree::Plane& plane, Direction direction ); 00253 00254 virtual ErrorCode initialize( BSPTree* tool, 00255 EntityHandle root, 00256 const double* point = 0 ); 00257 00258 public: 00259 00260 BSPTreeIter() : treeTool(0), childVect(2) {} 00261 virtual ~BSPTreeIter() {} 00262 00263 00264 BSPTree* tool() const 00265 { return treeTool; } 00266 00268 EntityHandle handle() const 00269 { return mStack.back(); } 00270 00272 unsigned depth() const 00273 { return mStack.size(); } 00274 00279 virtual ErrorCode step( Direction direction ); 00280 00285 ErrorCode step() { return step(RIGHT); } 00286 00291 ErrorCode back() { return step(LEFT); } 00292 00293 00295 ErrorCode get_parent_split_plane( BSPTree::Plane& plane ) const; 00296 00298 virtual double volume() const; 00299 00315 virtual bool intersect_ray( const double ray_point[3], 00316 const double ray_vect[3], 00317 double& t_enter, double& t_exit ) const; 00318 00321 bool is_sibling( const BSPTreeIter& other_leaf ) const; 00322 00325 bool is_sibling( EntityHandle other_leaf ) const; 00326 00331 bool sibling_is_forward( ) const; 00332 00334 virtual ErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const; 00335 }; 00336 00340 class BSPTreeBoxIter : public BSPTreeIter 00341 { 00342 private: 00343 00344 double leafCoords[8][3]; 00345 struct Corners { double coords[4][3]; }; 00346 std::vector<Corners> stackData; 00347 00348 protected: 00349 00350 virtual ErrorCode step_to_first_leaf( Direction direction ); 00351 00352 virtual ErrorCode up(); 00353 virtual ErrorCode down( const BSPTree::Plane& plane, Direction direction ); 00354 00355 virtual ErrorCode initialize( BSPTree* tool, 00356 EntityHandle root, 00357 const double* point = 0 ); 00358 public: 00359 00360 00362 enum SideBits { B0154 = 0x33, 00363 B1265 = 0x66, 00364 B2376 = 0xCC, 00365 B3047 = 0x99, 00366 B3210 = 0x0F, 00367 B4567 = 0xF0 00368 }; 00369 00370 static SideBits side_above_plane( const double hex_coords[8][3], 00371 const BSPTree::Plane& plane ); 00372 00373 static SideBits side_on_plane( const double hex_coords[8][3], 00374 const BSPTree::Plane& plane ); 00375 00376 static SideBits opposite_face( const SideBits& bits ) 00377 { return (SideBits)((~bits) & 0xFF); } 00378 00379 static ErrorCode face_corners( const SideBits face, 00380 const double hex_corners[8][3], 00381 double face_corners_out[4][3] ); 00382 00387 virtual ErrorCode step( Direction direction ); 00388 00393 ErrorCode step() { return BSPTreeIter::step(); } 00394 00399 ErrorCode back() { return BSPTreeIter::back(); } 00400 00402 ErrorCode get_box_corners( double coords[8][3] ) const; 00403 00405 double volume() const; 00406 00408 enum XSect { MISS = 0, SPLIT = 1, NONHEX = -1 }; 00409 XSect splits( const BSPTree::Plane& plane ) const; 00410 00412 bool intersects( const BSPTree::Plane& plane ) const; 00413 00429 bool intersect_ray( const double ray_point[3], 00430 const double ray_vect[3], 00431 double& t_enter, double& t_exit ) const; 00432 00441 ErrorCode sibling_side( SideBits& side_out ) const; 00442 00456 ErrorCode get_neighbors( SideBits side, 00457 std::vector<BSPTreeBoxIter>& results, 00458 double epsilon = 0.0 ) const; 00459 00461 ErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const; 00462 }; 00463 00464 } // namespace moab 00465 00466 #endif