MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "moab/GeomQueryTool.hpp" 00002 00003 #ifdef WIN32 /* windows */ 00004 #define _USE_MATH_DEFINES // For M_PI 00005 #endif 00006 #include <string> 00007 #include <iostream> 00008 #include <fstream> 00009 #include <sstream> 00010 #include <limits> 00011 #include <algorithm> 00012 #include <set> 00013 00014 #include <ctype.h> 00015 #include <string.h> 00016 #include <stdlib.h> 00017 #include <stdio.h> 00018 00019 #include "moab/OrientedBoxTreeTool.hpp" 00020 00021 const bool debug = false; 00022 #ifdef __DEBUG 00023 debug = true; 00024 #endif 00025 00026 namespace moab 00027 { 00028 00029 /** \class FindVolume_IntRegCtxt 00030 * 00031 * 00032 *\brief An intersection context used for finding a volume 00033 * 00034 * This context is used to find the nearest intersection location 00035 * and is intended for use with a global surface tree from 00036 * the GeomTopoTool. 00037 * 00038 * The behavior of this context is relatively simple in that it 00039 * returns only one intersection distance, surface, and facet. 00040 * Intersections of any orientation are accepted. The positive 00041 * value of the search window is limited to the current nearest 00042 * intersection distance. 00043 * 00044 */ 00045 00046 class FindVolumeIntRegCtxt : public OrientedBoxTreeTool::IntRegCtxt 00047 { 00048 00049 public: 00050 // Constructor 00051 FindVolumeIntRegCtxt() 00052 { 00053 // initialize return vectors 00054 // only one hit is returned in this context 00055 intersections.push_back( std::numeric_limits< double >::max() ); 00056 sets.push_back( 0 ); 00057 facets.push_back( 0 ); 00058 } 00059 00060 ErrorCode register_intersection( EntityHandle set, EntityHandle tri, double dist, 00061 OrientedBoxTreeTool::IntersectSearchWindow& search_win, 00062 GeomUtil::intersection_type /*it*/ ) 00063 { 00064 // update dist, set, and triangle hit if 00065 // we found a new minimum distance 00066 double abs_dist = fabs( dist ); 00067 if( abs_dist < fabs( intersections[0] ) ) 00068 { 00069 intersections[0] = dist; 00070 sets[0] = set; 00071 facets[0] = tri; 00072 00073 // narrow search window based on the hit distance 00074 pos = abs_dist; 00075 neg = -abs_dist; 00076 search_win.first = &pos; 00077 search_win.second = &neg; 00078 } 00079 00080 return MB_SUCCESS; 00081 } 00082 00083 // storage for updated window values during search 00084 double pos; 00085 double neg; 00086 }; 00087 00088 /** \class GQT_IntRegCtxt 00089 * 00090 *\brief An implementation of an Intersection Registration Context for use GQT ray-firing 00091 * 00092 * This context uses a variety of tests and conditions to confirm whether or 00093 * not to accumulate an intersection, to ensure robustness for ray firing. 00094 * 00095 * This context only accumulates intersections that are oriented parallel to 00096 * the 'desiredOrient', if provided, with respect to 'geomVol', using 00097 * information in the in 'senseTag'. 00098 * 00099 * This context only accumulates a single intersection out of a set of 00100 * multiple intersections that fall in the same 'neighborhood', where a 00101 * 'neighborhood' is defined as facets that share edges or vertices. 00102 * 00103 * This context only accumulates piercing intersections. This is relevant 00104 * for intersections that are found to be on an edge or vertex by the 00105 * Plucker test. Such intersections are piercing if the ray has the same 00106 * orientation w.r.t. to all fecets that share that edge or vertex. 00107 * 00108 * This context tests intersections against a list of 'prevFacets' to 00109 * prevent a ray from crossing the same facet more than once. The user is 00110 * responsible for ensuring that this list is reset when appropriate. 00111 * 00112 * This context accumulates all intersections within 'tol' of the 00113 * start of the ray and if the number of intersections within the 00114 * 'tol' of the ray start point is less than 'minTolInt', the next 00115 * closest intersection. If the desired result is only the closest 00116 * intersection, 'minTolInt' should be 0. This function will return all 00117 * intersections, regardless of distance from the start of the ray, if 00118 * 'minTolInt' is negative. 00119 * 00120 */ 00121 00122 class GQT_IntRegCtxt : public OrientedBoxTreeTool::IntRegCtxt 00123 { 00124 00125 private: 00126 // Input 00127 OrientedBoxTreeTool* tool; 00128 const CartVect ray_origin; 00129 const CartVect ray_direction; 00130 const double tol; /* used for box.intersect_ray, radius of 00131 neighborhood for adjacent triangles, 00132 and old mode of add_intersection */ 00133 const int minTolInt; 00134 00135 // Optional Input - to screen RTIs by orientation and edge/node intersection 00136 const EntityHandle* rootSet; /* used for sphere_intersect */ 00137 const EntityHandle* geomVol; /* used for determining surface sense */ 00138 const Tag* senseTag; /* allows screening by triangle orientation. 00139 both geomVol and senseTag must be used together. */ 00140 const int* desiredOrient; /* points to desired orientation of ray with 00141 respect to surf normal, if this feature is used. 00142 Must point to -1 (reverse) or 1 (forward). 00143 geomVol and senseTag are needed for this feature */ 00144 00145 // Optional Input - to avoid returning these as RTIs 00146 const std::vector< EntityHandle >* prevFacets; /* intersections on these triangles 00147 will not be returned */ 00148 00149 // Other Variables 00150 std::vector< std::vector< EntityHandle > > neighborhoods; 00151 std::vector< EntityHandle > neighborhood; 00152 00153 void add_intersection( EntityHandle set, EntityHandle tri, double dist, 00154 OrientedBoxTreeTool::IntersectSearchWindow& search_win ); 00155 void append_intersection( EntityHandle set, EntityHandle facet, double dist ); 00156 void set_intersection( int len_idx, EntityHandle set, EntityHandle facet, double dist ); 00157 void add_mode1_intersection( EntityHandle set, EntityHandle facet, double dist, 00158 OrientedBoxTreeTool::IntersectSearchWindow& search_win ); 00159 bool edge_node_piercing_intersect( const EntityHandle tri, const CartVect& ray_direction, 00160 const GeomUtil::intersection_type int_type, 00161 const std::vector< EntityHandle >& close_tris, 00162 const std::vector< int >& close_senses, const Interface* MBI, 00163 std::vector< EntityHandle >* neighborhood_tris = 0 ); 00164 00165 bool in_prevFacets( const EntityHandle tri ); 00166 bool in_neighborhoods( const EntityHandle tri ); 00167 00168 public: 00169 GQT_IntRegCtxt( OrientedBoxTreeTool* obbtool, const double ray_point[3], const double ray_dir[3], double tolerance, 00170 int min_tolerance_intersections, const EntityHandle* root_set, const EntityHandle* geom_volume, 00171 const Tag* sense_tag, const int* desired_orient, const std::vector< EntityHandle >* prev_facets ) 00172 : tool( obbtool ), ray_origin( ray_point ), ray_direction( ray_dir ), tol( tolerance ), 00173 minTolInt( min_tolerance_intersections ), rootSet( root_set ), geomVol( geom_volume ), senseTag( sense_tag ), 00174 desiredOrient( desired_orient ), prevFacets( prev_facets ){ 00175 00176 }; 00177 00178 virtual ErrorCode register_intersection( EntityHandle set, EntityHandle triangle, double distance, 00179 OrientedBoxTreeTool::IntersectSearchWindow&, 00180 GeomUtil::intersection_type int_type ); 00181 00182 virtual ErrorCode update_orient( EntityHandle set, int* surfTriOrient ); 00183 virtual const int* getDesiredOrient() 00184 { 00185 return desiredOrient; 00186 }; 00187 }; 00188 00189 ErrorCode GQT_IntRegCtxt::update_orient( EntityHandle set, int* surfTriOrient ) 00190 { 00191 00192 ErrorCode rval; 00193 00194 // Get desired orientation of surface wrt volume. Use this to return only 00195 // exit or entrance intersections. 00196 if( geomVol && senseTag && desiredOrient && surfTriOrient ) 00197 { 00198 if( 1 != *desiredOrient && -1 != *desiredOrient ) 00199 { std::cerr << "error: desired orientation must be 1 (forward) or -1 (reverse)" << std::endl; } 00200 EntityHandle vols[2]; 00201 rval = tool->get_moab_instance()->tag_get_data( *senseTag, &set, 1, vols ); 00202 assert( MB_SUCCESS == rval ); 00203 if( MB_SUCCESS != rval ) return rval; 00204 if( vols[0] == vols[1] ) 00205 { 00206 std::cerr << "error: surface has positive and negative sense wrt same volume" << std::endl; 00207 return MB_FAILURE; 00208 } 00209 // surfTriOrient will be used by plucker_ray_tri_intersect to avoid 00210 // intersections with wrong orientation. 00211 if( *geomVol == vols[0] ) { *surfTriOrient = *desiredOrient * 1; } 00212 else if( *geomVol == vols[1] ) 00213 { 00214 *surfTriOrient = *desiredOrient * ( -1 ); 00215 } 00216 else 00217 { 00218 assert( false ); 00219 return MB_FAILURE; 00220 } 00221 } 00222 00223 return MB_SUCCESS; 00224 } 00225 00226 bool GQT_IntRegCtxt::in_prevFacets( const EntityHandle tri ) 00227 { 00228 return ( prevFacets && ( ( *prevFacets ).end() != find( ( *prevFacets ).begin(), ( *prevFacets ).end(), tri ) ) ); 00229 } 00230 00231 bool GQT_IntRegCtxt::in_neighborhoods( const EntityHandle tri ) 00232 { 00233 bool same_neighborhood = false; 00234 for( unsigned i = 0; i < neighborhoods.size(); ++i ) 00235 { 00236 if( neighborhoods[i].end() != find( neighborhoods[i].begin(), neighborhoods[i].end(), tri ) ) 00237 { 00238 same_neighborhood = true; 00239 continue; 00240 } 00241 } 00242 return same_neighborhood; 00243 } 00244 00245 /**\brief Determine if a ray-edge/node intersection is glancing or piercing. 00246 * This function avoids asking for upward adjacencies to prevent their 00247 * creation. 00248 *\param tri The intersected triangle 00249 *\param ray_dir The direction of the ray 00250 *\param int_type The type of intersection (EDGE0, EDGE1, NODE2, ...) 00251 *\param close_tris Vector of triangles in the proximity of the intersection 00252 *\param close_senses Vector of surface senses for tris in the proximity of 00253 * the intersection 00254 *\param neighborhood Vector of triangles in the topological neighborhood of the intersection 00255 *\return True if piercing, false otherwise. 00256 */ 00257 bool GQT_IntRegCtxt::edge_node_piercing_intersect( const EntityHandle tri, const CartVect& ray_dir, 00258 const GeomUtil::intersection_type int_type, 00259 const std::vector< EntityHandle >& close_tris, 00260 const std::vector< int >& close_senses, const Interface* MBI, 00261 std::vector< EntityHandle >* neighborhood_tris ) 00262 { 00263 00264 // get the node of the triangle 00265 const EntityHandle* conn = NULL; 00266 int len = 0; 00267 ErrorCode rval = MBI->get_connectivity( tri, conn, len ); 00268 if( MB_SUCCESS != rval || 3 != len ) return MB_FAILURE; 00269 00270 // get adjacent tris (and keep their corresponding senses) 00271 std::vector< EntityHandle > adj_tris; 00272 std::vector< int > adj_senses; 00273 00274 // node intersection 00275 if( GeomUtil::NODE0 == int_type || GeomUtil::NODE1 == int_type || GeomUtil::NODE2 == int_type ) 00276 { 00277 00278 // get the intersected node 00279 EntityHandle node; 00280 if( GeomUtil::NODE0 == int_type ) 00281 node = conn[0]; 00282 else if( GeomUtil::NODE1 == int_type ) 00283 node = conn[1]; 00284 else 00285 node = conn[2]; 00286 00287 // get tris adjacent to node 00288 for( unsigned i = 0; i < close_tris.size(); ++i ) 00289 { 00290 const EntityHandle* con = NULL; 00291 rval = MBI->get_connectivity( close_tris[i], con, len ); 00292 if( MB_SUCCESS != rval || 3 != len ) return MB_FAILURE; 00293 00294 if( node == con[0] || node == con[1] || node == con[2] ) 00295 { 00296 adj_tris.push_back( close_tris[i] ); 00297 adj_senses.push_back( close_senses[i] ); 00298 } 00299 } 00300 if( adj_tris.empty() ) 00301 { 00302 std::cerr << "error: no tris are adjacent to the node" << std::endl; 00303 return MB_FAILURE; 00304 } 00305 // edge intersection 00306 } 00307 else if( GeomUtil::EDGE0 == int_type || GeomUtil::EDGE1 == int_type || GeomUtil::EDGE2 == int_type ) 00308 { 00309 00310 // get the endpoints of the edge 00311 EntityHandle endpts[2]; 00312 if( GeomUtil::EDGE0 == int_type ) 00313 { 00314 endpts[0] = conn[0]; 00315 endpts[1] = conn[1]; 00316 } 00317 else if( GeomUtil::EDGE1 == int_type ) 00318 { 00319 endpts[0] = conn[1]; 00320 endpts[1] = conn[2]; 00321 } 00322 else 00323 { 00324 endpts[0] = conn[2]; 00325 endpts[1] = conn[0]; 00326 } 00327 00328 // get tris adjacent to edge 00329 for( unsigned i = 0; i < close_tris.size(); ++i ) 00330 { 00331 const EntityHandle* con = NULL; 00332 rval = MBI->get_connectivity( close_tris[i], con, len ); 00333 if( MB_SUCCESS != rval || 3 != len ) return MB_FAILURE; 00334 00335 // check both orientations in case close_tris are not on the same surface 00336 if( ( endpts[0] == con[0] && endpts[1] == con[1] ) || ( endpts[0] == con[1] && endpts[1] == con[0] ) || 00337 ( endpts[0] == con[1] && endpts[1] == con[2] ) || ( endpts[0] == con[2] && endpts[1] == con[1] ) || 00338 ( endpts[0] == con[2] && endpts[1] == con[0] ) || ( endpts[0] == con[0] && endpts[1] == con[2] ) ) 00339 { 00340 adj_tris.push_back( close_tris[i] ); 00341 adj_senses.push_back( close_senses[i] ); 00342 } 00343 } 00344 // In a 2-manifold each edge is adjacent to exactly 2 tris 00345 if( 2 != adj_tris.size() ) 00346 { 00347 std::cerr << "error: edge of a manifold must be topologically adjacent to exactly 2 tris" << std::endl; 00348 MBI->list_entities( endpts, 2 ); 00349 return true; 00350 } 00351 } 00352 else 00353 { 00354 std::cerr << "error: special case not an node/edge intersection" << std::endl; 00355 return MB_FAILURE; 00356 } 00357 00358 // The close tris were in proximity to the intersection. The adj_tris are 00359 // topologically adjacent to the intersection (the neighborhood). 00360 if( neighborhood_tris ) ( *neighborhood_tris ).assign( adj_tris.begin(), adj_tris.end() ); 00361 00362 // determine glancing/piercing 00363 // If a desired_orientation was used in this call to ray_intersect_sets, 00364 // the plucker_ray_tri_intersect will have already used it. For a piercing 00365 // intersection, the normal of all tris must have the same orientation. 00366 int sign = 0; 00367 for( unsigned i = 0; i < adj_tris.size(); ++i ) 00368 { 00369 const EntityHandle* con = NULL; 00370 rval = MBI->get_connectivity( adj_tris[i], con, len ); 00371 if( MB_SUCCESS != rval || 3 != len ) return MB_FAILURE; 00372 CartVect coords[3]; 00373 rval = MBI->get_coords( con, len, coords[0].array() ); 00374 if( MB_SUCCESS != rval ) return MB_FAILURE; 00375 00376 // get normal of triangle 00377 CartVect v0 = coords[1] - coords[0]; 00378 CartVect v1 = coords[2] - coords[0]; 00379 CartVect norm = adj_senses[i] * ( v0 * v1 ); 00380 double dot_prod = norm % ray_dir; 00381 00382 // if the sign has not yet been decided, choose it 00383 if( 0 == sign && 0 != dot_prod ) 00384 { 00385 if( 0 < dot_prod ) 00386 sign = 1; 00387 else 00388 sign = -1; 00389 } 00390 00391 // intersection is glancing if tri and ray do not point in same direction 00392 // for every triangle 00393 if( 0 != sign && 0 > sign * dot_prod ) return false; 00394 } 00395 return true; 00396 } 00397 00398 ErrorCode GQT_IntRegCtxt::register_intersection( EntityHandle set, EntityHandle t, double int_dist, 00399 OrientedBoxTreeTool::IntersectSearchWindow& search_win, 00400 GeomUtil::intersection_type int_type ) 00401 { 00402 ErrorCode rval; 00403 00404 // Do not accept intersections if they are in the vector of previously intersected 00405 // facets. 00406 if( in_prevFacets( t ) ) return MB_SUCCESS; 00407 00408 // Do not accept intersections if they are in the neighborhood of previous 00409 // intersections. 00410 if( in_neighborhoods( t ) ) return MB_SUCCESS; 00411 00412 neighborhood.clear(); 00413 00414 // Handle special case of edge/node intersection. Accept piercing 00415 // intersections and reject glancing intersections. 00416 // The edge_node_intersection function needs to know surface sense wrt volume. 00417 // A less-robust implementation could work without sense information. 00418 // Would it ever be useful to accept a glancing intersection? 00419 if( GeomUtil::INTERIOR != int_type && rootSet && geomVol && senseTag ) 00420 { 00421 // get triangles in the proximity of the intersection 00422 std::vector< EntityHandle > close_tris; 00423 std::vector< int > close_senses; 00424 rval = tool->get_close_tris( ray_origin + int_dist * ray_direction, tol, rootSet, geomVol, senseTag, close_tris, 00425 close_senses ); 00426 00427 if( MB_SUCCESS != rval ) return rval; 00428 00429 if( !edge_node_piercing_intersect( t, ray_direction, int_type, close_tris, close_senses, 00430 tool->get_moab_instance(), &neighborhood ) ) 00431 return MB_SUCCESS; 00432 } 00433 else 00434 { 00435 neighborhood.push_back( t ); 00436 } 00437 00438 // NOTE: add_intersection may modify the 'neg_ray_len' and 'nonneg_ray_len' 00439 // members, which will affect subsequent calls to ray_tri_intersect 00440 // in this loop. 00441 add_intersection( set, t, int_dist, search_win ); 00442 00443 return MB_SUCCESS; 00444 } 00445 00446 void GQT_IntRegCtxt::append_intersection( EntityHandle set, EntityHandle facet, double dist ) 00447 { 00448 intersections.push_back( dist ); 00449 sets.push_back( set ); 00450 facets.push_back( facet ); 00451 neighborhoods.push_back( neighborhood ); 00452 return; 00453 } 00454 00455 void GQT_IntRegCtxt::set_intersection( int len_idx, EntityHandle set, EntityHandle facet, double dist ) 00456 { 00457 intersections[len_idx] = dist; 00458 sets[len_idx] = set; 00459 facets[len_idx] = facet; 00460 return; 00461 } 00462 00463 /* Mode 1: Used if neg_ray_len and nonneg_ray_len are specified 00464 variables used: nonneg_ray_len, neg_ray_len 00465 variables not used: min_tol_int, tol 00466 1) keep the closest nonneg intersection and one negative intersection, if closer 00467 */ 00468 void GQT_IntRegCtxt::add_mode1_intersection( EntityHandle set, EntityHandle facet, double dist, 00469 OrientedBoxTreeTool::IntersectSearchWindow& search_win ) 00470 { 00471 if( 2 != intersections.size() ) 00472 { 00473 intersections.resize( 2, 0 ); 00474 sets.resize( 2, 0 ); 00475 facets.resize( 2, 0 ); 00476 // must initialize this for comparison below 00477 intersections[0] = -std::numeric_limits< double >::max(); 00478 } 00479 00480 // negative case 00481 if( 0.0 > dist ) 00482 { 00483 set_intersection( 0, set, facet, dist ); 00484 search_win.second = &intersections[0]; 00485 // nonnegative case 00486 } 00487 else 00488 { 00489 set_intersection( 1, set, facet, dist ); 00490 search_win.first = &intersections[1]; 00491 // if the intersection is closer than the negative one, remove the negative one 00492 if( dist < -*( search_win.second ) ) 00493 { 00494 set_intersection( 0, 0, 0, -intersections[1] ); 00495 search_win.second = &intersections[0]; 00496 } 00497 } 00498 // std::cout << "add_intersection: dist = " << dist << " search_win.second=" << 00499 // *search_win.second 00500 // << " search_win.first=" << *search_win.first << std::endl; 00501 return; 00502 } 00503 00504 void GQT_IntRegCtxt::add_intersection( EntityHandle set, EntityHandle facet, double dist, 00505 OrientedBoxTreeTool::IntersectSearchWindow& search_win ) 00506 { 00507 00508 // Mode 1, detected by non-null neg_ray_len pointer 00509 // keep the closest nonneg intersection and one negative intersection, if closer 00510 if( search_win.second && search_win.first ) { return add_mode1_intersection( set, facet, dist, search_win ); } 00511 00512 // --------------------------------------------------------------------------- 00513 /* Mode 2: Used if neg_ray_len is not specified 00514 variables used: min_tol_int, tol, search_win.first 00515 variables not used: neg_ray_len 00516 1) if(min_tol_int<0) return all intersections 00517 2) otherwise return all inside tolerance and unless there are >min_tol_int 00518 inside of tolerance, return the closest outside of tolerance */ 00519 // Mode 2 00520 // If minTolInt is less than zero, return all intersections 00521 if( minTolInt < 0 && dist > -tol ) 00522 { 00523 append_intersection( set, facet, dist ); 00524 neighborhoods.push_back( neighborhood ); 00525 return; 00526 } 00527 00528 // Check if the 'len' pointer is pointing into the intersection 00529 // list. If this is the case, then the list contains, at that 00530 // location, an intersection greater than the tolerance away from 00531 // the base point of the ray. 00532 int len_idx = -1; 00533 if( search_win.first && search_win.first >= &intersections[0] && 00534 search_win.first < &intersections[0] + intersections.size() ) 00535 len_idx = search_win.first - &intersections[0]; 00536 00537 // If the intersection is within tol of the ray base point, we 00538 // always add it to the list. 00539 if( dist <= tol ) 00540 { 00541 // If the list contains an intersection outside the tolerance... 00542 if( len_idx >= 0 ) 00543 { 00544 // If we no longer want an intersection outside the tolerance, 00545 // remove it. 00546 if( (int)intersections.size() >= minTolInt ) 00547 { 00548 set_intersection( len_idx, set, facet, dist ); 00549 // From now on, we want only intersections within the tolerance, 00550 // so update length accordingly 00551 search_win.first = &tol; 00552 } 00553 // Otherwise appended to the list and update pointer 00554 else 00555 { 00556 append_intersection( set, facet, dist ); 00557 search_win.first = &intersections[len_idx]; 00558 } 00559 } 00560 // Otherwise just append it 00561 else 00562 { 00563 append_intersection( set, facet, dist ); 00564 // If we have all the intersections we want, set 00565 // length such that we will only find further intersections 00566 // within the tolerance 00567 if( (int)intersections.size() >= minTolInt ) search_win.first = &tol; 00568 } 00569 } 00570 // Otherwise the intersection is outside the tolerance 00571 // If we already have an intersection outside the tolerance and 00572 // this one is closer, replace the existing one with this one. 00573 else if( len_idx >= 0 ) 00574 { 00575 if( dist <= *search_win.first ) { set_intersection( len_idx, set, facet, dist ); } 00576 } 00577 // Otherwise if we want an intersection outside the tolerance 00578 // and don'thave one yet, add it. 00579 else if( (int)intersections.size() < minTolInt ) 00580 { 00581 append_intersection( set, facet, dist ); 00582 // update length. this is currently the closest intersection, so 00583 // only want further intersections that are closer than this one. 00584 search_win.first = &intersections.back(); 00585 } 00586 } 00587 00588 GeomQueryTool::GeomQueryTool( Interface* impl, bool find_geomsets, EntityHandle modelRootSet, bool p_rootSets_vector, 00589 bool restore_rootSets, bool trace_counting, double overlap_thickness, 00590 double numerical_precision ) 00591 : owns_gtt( true ) 00592 { 00593 geomTopoTool = new GeomTopoTool( impl, find_geomsets, modelRootSet, p_rootSets_vector, restore_rootSets ); 00594 00595 senseTag = geomTopoTool->get_sense_tag(); 00596 00597 obbTreeTool = geomTopoTool->obb_tree(); 00598 MBI = geomTopoTool->get_moab_instance(); 00599 00600 counting = trace_counting; 00601 overlapThickness = overlap_thickness; 00602 numericalPrecision = numerical_precision; 00603 00604 // reset query counters 00605 n_pt_in_vol_calls = 0; 00606 n_ray_fire_calls = 0; 00607 } 00608 00609 GeomQueryTool::GeomQueryTool( GeomTopoTool* geomtopotool, bool trace_counting, double overlap_thickness, 00610 double numerical_precision ) 00611 : owns_gtt( false ) 00612 { 00613 00614 geomTopoTool = geomtopotool; 00615 00616 senseTag = geomTopoTool->get_sense_tag(); 00617 00618 obbTreeTool = geomTopoTool->obb_tree(); 00619 MBI = geomTopoTool->get_moab_instance(); 00620 00621 counting = trace_counting; 00622 overlapThickness = overlap_thickness; 00623 numericalPrecision = numerical_precision; 00624 00625 // reset query counters 00626 n_pt_in_vol_calls = 0; 00627 n_ray_fire_calls = 0; 00628 } 00629 00630 GeomQueryTool::~GeomQueryTool() 00631 { 00632 if( owns_gtt ) { delete geomTopoTool; } 00633 } 00634 00635 ErrorCode GeomQueryTool::initialize() 00636 { 00637 00638 ErrorCode rval; 00639 00640 rval = geomTopoTool->find_geomsets();MB_CHK_SET_ERR( rval, "Failed to find geometry sets" ); 00641 00642 rval = geomTopoTool->setup_implicit_complement();MB_CHK_SET_ERR( rval, "Couldn't setup the implicit complement" ); 00643 00644 rval = geomTopoTool->construct_obb_trees();MB_CHK_SET_ERR( rval, "Failed to construct OBB trees" ); 00645 00646 return MB_SUCCESS; 00647 } 00648 00649 void GeomQueryTool::RayHistory::reset() 00650 { 00651 prev_facets.clear(); 00652 } 00653 00654 void GeomQueryTool::RayHistory::reset_to_last_intersection() 00655 { 00656 00657 if( prev_facets.size() > 1 ) 00658 { 00659 prev_facets[0] = prev_facets.back(); 00660 prev_facets.resize( 1 ); 00661 } 00662 } 00663 00664 void GeomQueryTool::RayHistory::rollback_last_intersection() 00665 { 00666 if( prev_facets.size() ) prev_facets.pop_back(); 00667 } 00668 00669 ErrorCode GeomQueryTool::RayHistory::get_last_intersection( EntityHandle& last_facet_hit ) const 00670 { 00671 if( prev_facets.size() > 0 ) 00672 { 00673 last_facet_hit = prev_facets.back(); 00674 return MB_SUCCESS; 00675 } 00676 else 00677 { 00678 return MB_ENTITY_NOT_FOUND; 00679 } 00680 } 00681 00682 bool GeomQueryTool::RayHistory::in_history( EntityHandle ent ) const 00683 { 00684 return std::find( prev_facets.begin(), prev_facets.end(), ent ) != prev_facets.end(); 00685 } 00686 00687 void GeomQueryTool::RayHistory::add_entity( EntityHandle ent ) 00688 { 00689 prev_facets.push_back( ent ); 00690 } 00691 00692 ErrorCode GeomQueryTool::ray_fire( const EntityHandle volume, const double point[3], const double dir[3], 00693 EntityHandle& next_surf, double& next_surf_dist, RayHistory* history, 00694 double user_dist_limit, int ray_orientation, OrientedBoxTreeTool::TrvStats* stats ) 00695 { 00696 00697 // take some stats that are independent of nps 00698 if( counting ) 00699 { 00700 ++n_ray_fire_calls; 00701 if( 0 == n_ray_fire_calls % 10000000 ) 00702 { std::cout << "n_ray_fires=" << n_ray_fire_calls << " n_pt_in_vols=" << n_pt_in_vol_calls << std::endl; } 00703 } 00704 00705 if( debug ) 00706 { 00707 std::cout << "ray_fire:" 00708 << " xyz=" << point[0] << " " << point[1] << " " << point[2] << " uvw=" << dir[0] << " " << dir[1] 00709 << " " << dir[2] << " entity_handle=" << volume << std::endl; 00710 } 00711 00712 const double huge_val = std::numeric_limits< double >::max(); 00713 double dist_limit = huge_val; 00714 if( user_dist_limit > 0 ) dist_limit = user_dist_limit; 00715 00716 // don't recreate these every call 00717 std::vector< double > dists; 00718 std::vector< EntityHandle > surfs; 00719 std::vector< EntityHandle > facets; 00720 00721 EntityHandle root; 00722 ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the obb tree root of the volume" ); 00723 00724 // check behind the ray origin for intersections 00725 double neg_ray_len; 00726 if( 0 == overlapThickness ) { neg_ray_len = -numericalPrecision; } 00727 else 00728 { 00729 neg_ray_len = -overlapThickness; 00730 } 00731 00732 // optionally, limit the nonneg_ray_len with the distance to next collision. 00733 double nonneg_ray_len = dist_limit; 00734 00735 // the nonneg_ray_len should not be less than -neg_ray_len, or an overlap 00736 // may be missed due to optimization within ray_intersect_sets 00737 if( nonneg_ray_len < -neg_ray_len ) nonneg_ray_len = -neg_ray_len; 00738 if( 0 > nonneg_ray_len || 0 <= neg_ray_len ) { MB_SET_ERR( MB_FAILURE, "Incorrect ray length provided" ); } 00739 00740 // min_tolerance_intersections is passed but not used in this call 00741 const int min_tolerance_intersections = 0; 00742 00743 // numericalPrecision is used for box.intersect_ray and find triangles in the 00744 // neighborhood of edge/node intersections. 00745 GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), point, dir, numericalPrecision, min_tolerance_intersections, 00746 &root, &volume, &senseTag, &ray_orientation, 00747 history ? &( history->prev_facets ) : NULL ); 00748 00749 OrientedBoxTreeTool::IntersectSearchWindow search_win( &nonneg_ray_len, &neg_ray_len ); 00750 rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, point, dir, 00751 search_win, int_reg_ctxt, stats ); 00752 00753 MB_CHK_SET_ERR( rval, "Ray query failed" ); 00754 00755 // If no distances are returned, the particle is lost unless the physics limit 00756 // is being used. If the physics limit is being used, there is no way to tell 00757 // if the particle is lost. To avoid ambiguity, DO NOT use the distance limit 00758 // unless you know lost particles do not occur. 00759 if( dists.empty() ) 00760 { 00761 next_surf = 0; 00762 if( debug ) { std::cout << " next_surf=0 dist=(undef)" << std::endl; } 00763 return MB_SUCCESS; 00764 } 00765 00766 // Assume that a (neg, nonneg) pair of RTIs could be returned, 00767 // however, only one or the other may exist. dists[] may be populated, but 00768 // intersections are ONLY indicated by nonzero surfs[] and facets[]. 00769 if( 2 != dists.size() || 2 != facets.size() ) { MB_SET_ERR( MB_FAILURE, "Incorrect number of facets/distances" ); } 00770 if( 0.0 < dists[0] || 0.0 > dists[1] ) { MB_SET_ERR( MB_FAILURE, "Invalid intersection distance signs" ); } 00771 00772 // If both negative and nonnegative RTIs are returned, the negative RTI must 00773 // closer to the origin. 00774 if( ( 0 != facets[0] && 0 != facets[1] ) && ( -dists[0] > dists[1] ) ) 00775 { MB_SET_ERR( MB_FAILURE, "Invalid intersection distance values" ); } 00776 00777 // If an RTI is found at negative distance, perform a PMT to see if the 00778 // particle is inside an overlap. 00779 int exit_idx = -1; 00780 if( 0 != facets[0] ) 00781 { 00782 // get the next volume 00783 std::vector< EntityHandle > vols; 00784 EntityHandle nx_vol; 00785 rval = MBI->get_parent_meshsets( surfs[0], vols );MB_CHK_SET_ERR( rval, "Failed to get the parent meshsets" ); 00786 if( 2 != vols.size() ) { MB_SET_ERR( MB_FAILURE, "Invaid number of parent volumes found" ); } 00787 if( vols.front() == volume ) { nx_vol = vols.back(); } 00788 else 00789 { 00790 nx_vol = vols.front(); 00791 } 00792 // Check to see if the point is actually in the next volume. 00793 // The list of previous facets is used to topologically identify the 00794 // "on_boundary" result of the PMT. This avoids a test that uses proximity 00795 // (a tolerance). 00796 int result; 00797 rval = point_in_volume( nx_vol, point, result, dir, history );MB_CHK_SET_ERR( rval, "Point in volume query failed" ); 00798 if( 1 == result ) exit_idx = 0; 00799 } 00800 00801 // if the negative distance is not the exit, try the nonnegative distance 00802 if( -1 == exit_idx && 0 != facets[1] ) exit_idx = 1; 00803 00804 // if the exit index is still unknown, the particle is lost 00805 if( -1 == exit_idx ) 00806 { 00807 next_surf = 0; 00808 if( debug ) { std::cout << "next surf hit = 0, dist = (undef)" << std::endl; } 00809 return MB_SUCCESS; 00810 } 00811 00812 // return the intersection 00813 next_surf = surfs[exit_idx]; 00814 next_surf_dist = ( 0 > dists[exit_idx] ? 0 : dists[exit_idx] ); 00815 00816 if( history ) { history->prev_facets.push_back( facets[exit_idx] ); } 00817 00818 if( debug ) 00819 { 00820 if( 0 > dists[exit_idx] ) { std::cout << " OVERLAP track length=" << dists[exit_idx] << std::endl; } 00821 std::cout << " next_surf = " << next_surf // todo: use geomtopotool to get id by entity handle 00822 << ", dist = " << next_surf_dist << " new_pt="; 00823 for( int i = 0; i < 3; ++i ) 00824 { 00825 std::cout << point[i] + dir[i] * next_surf_dist << " "; 00826 } 00827 std::cout << std::endl; 00828 } 00829 00830 return MB_SUCCESS; 00831 } 00832 00833 ErrorCode GeomQueryTool::point_in_volume( const EntityHandle volume, const double xyz[3], int& result, 00834 const double* uvw, const RayHistory* history ) 00835 { 00836 // take some stats that are independent of nps 00837 if( counting ) ++n_pt_in_vol_calls; 00838 00839 // early fail for piv - see if point inside the root level obb 00840 // if its not even in the box dont bother doing anything else 00841 ErrorCode rval = point_in_box( volume, xyz, result ); 00842 if( !result ) 00843 { 00844 result = 0; 00845 return MB_SUCCESS; 00846 } 00847 00848 // get OBB Tree for volume 00849 EntityHandle root; 00850 rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to find the volume's obb tree root" ); 00851 00852 // Don't recreate these every call. These cannot be the same as the ray_fire 00853 // vectors because both are used simultaneously. 00854 std::vector< double > dists; 00855 std::vector< EntityHandle > surfs; 00856 std::vector< EntityHandle > facets; 00857 std::vector< int > dirs; 00858 00859 // if uvw is not given or is full of zeros, use a random direction 00860 double u = 0, v = 0, w = 0; 00861 00862 if( uvw ) 00863 { 00864 u = uvw[0]; 00865 v = uvw[1], w = uvw[2]; 00866 } 00867 00868 if( u == 0 && v == 0 && w == 0 ) 00869 { 00870 u = rand(); 00871 v = rand(); 00872 w = rand(); 00873 const double magnitude = sqrt( u * u + v * v + w * w ); 00874 u /= magnitude; 00875 v /= magnitude; 00876 w /= magnitude; 00877 } 00878 00879 const double ray_direction[] = { u, v, w }; 00880 00881 // if overlaps, ray must be cast to infinity and all RTIs must be returned 00882 const double large = 1e15; 00883 double ray_length = large; 00884 00885 // If overlaps occur, the pt is inside if traveling along the ray from the 00886 // origin, there are ever more exits than entrances. In lieu of implementing 00887 // that, all intersections to infinity are required if overlaps occur (expensive) 00888 int min_tolerance_intersections; 00889 if( 0 != overlapThickness ) 00890 { 00891 min_tolerance_intersections = -1; 00892 // only the first intersection is needed if overlaps do not occur (cheap) 00893 } 00894 else 00895 { 00896 min_tolerance_intersections = 1; 00897 } 00898 00899 // Get intersection(s) of forward and reverse orientation. Do not return 00900 // glancing intersections or previous facets. 00901 GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), xyz, ray_direction, numericalPrecision, 00902 min_tolerance_intersections, &root, &volume, &senseTag, NULL, 00903 history ? &( history->prev_facets ) : NULL ); 00904 00905 OrientedBoxTreeTool::IntersectSearchWindow search_win( &ray_length, (double*)NULL ); 00906 rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, xyz, 00907 ray_direction, search_win, int_reg_ctxt );MB_CHK_SET_ERR( rval, "Ray fire query failed" ); 00908 00909 // determine orientation of all intersections 00910 // 1 for entering, 0 for leaving, -1 for tangent 00911 // Tangent intersections are not returned from ray_tri_intersect. 00912 dirs.resize( dists.size() ); 00913 for( unsigned i = 0; i < dists.size(); ++i ) 00914 { 00915 rval = boundary_case( volume, dirs[i], u, v, w, facets[i], surfs[i] );MB_CHK_SET_ERR( rval, "Failed to resolve boundary case" ); 00916 } 00917 00918 // count all crossings 00919 if( 0 != overlapThickness ) 00920 { 00921 int sum = 0; 00922 for( unsigned i = 0; i < dirs.size(); ++i ) 00923 { 00924 if( 1 == dirs[i] ) 00925 sum += 1; // +1 for entering 00926 else if( 0 == dirs[i] ) 00927 sum -= 1; // -1 for leaving 00928 else if( -1 == dirs[i] ) 00929 { // 0 for tangent 00930 std::cout << "direction==tangent" << std::endl; 00931 sum += 0; 00932 } 00933 else 00934 { 00935 MB_SET_ERR( MB_FAILURE, "Error: unknown direction" ); 00936 } 00937 } 00938 00939 // inside/outside depends on the sum 00940 if( 0 < sum ) 00941 result = 0; // pt is outside (for all vols) 00942 else if( 0 > sum ) 00943 result = 1; // pt is inside (for all vols) 00944 else if( geomTopoTool->is_implicit_complement( volume ) ) 00945 result = 1; // pt is inside (for impl_compl_vol) 00946 else 00947 result = 0; // pt is outside (for all other vols) 00948 00949 // Only use the first crossing 00950 } 00951 else 00952 { 00953 if( dirs.empty() ) 00954 { 00955 result = 0; // pt is outside 00956 } 00957 else 00958 { 00959 int smallest = std::min_element( dists.begin(), dists.end() ) - dists.begin(); 00960 if( 1 == dirs[smallest] ) 00961 result = 0; // pt is outside 00962 else if( 0 == dirs[smallest] ) 00963 result = 1; // pt is inside 00964 else if( -1 == dirs[smallest] ) 00965 { 00966 // Should not be here because Plucker ray-triangle test does not 00967 // return coplanar rays as intersections. 00968 std::cout << "direction==tangent" << std::endl; 00969 result = -1; 00970 } 00971 else 00972 { 00973 MB_SET_ERR( MB_FAILURE, "Error: unknown direction" ); 00974 } 00975 } 00976 } 00977 00978 if( debug ) 00979 std::cout << "pt_in_vol: result=" << result << " xyz=" << xyz[0] << " " << xyz[1] << " " << xyz[2] 00980 << " uvw=" << u << " " << v << " " << w << " vol_id=" << volume 00981 << std::endl; // todo: use geomtopotool to get id by entity handle 00982 00983 return MB_SUCCESS; 00984 } 00985 00986 /** 00987 * \brief For the volume pointed to and the point wished to be tested, returns 00988 * whether the point is inside or outside the bounding box of the volume. 00989 * inside = 0, not inside, inside = 1, inside 00990 */ 00991 ErrorCode GeomQueryTool::point_in_box( EntityHandle volume, const double point[3], int& inside ) 00992 { 00993 double minpt[3]; 00994 double maxpt[3]; 00995 ErrorCode rval = geomTopoTool->get_bounding_coords( volume, minpt, maxpt );MB_CHK_SET_ERR( rval, "Failed to get the bounding coordinates of the volume" ); 00996 00997 // early exits 00998 if( point[0] > maxpt[0] || point[0] < minpt[0] ) 00999 { 01000 inside = 0; 01001 return rval; 01002 } 01003 if( point[1] > maxpt[1] || point[1] < minpt[1] ) 01004 { 01005 inside = 0; 01006 return rval; 01007 } 01008 if( point[2] > maxpt[2] || point[2] < minpt[2] ) 01009 { 01010 inside = 0; 01011 return rval; 01012 } 01013 inside = 1; 01014 return rval; 01015 } 01016 01017 ErrorCode GeomQueryTool::test_volume_boundary( const EntityHandle volume, const EntityHandle surface, 01018 const double xyz[3], const double uvw[3], int& result, 01019 const RayHistory* history ) 01020 { 01021 ErrorCode rval; 01022 int dir; 01023 01024 if( history && history->prev_facets.size() ) 01025 { 01026 // the current facet is already available 01027 rval = boundary_case( volume, dir, uvw[0], uvw[1], uvw[2], history->prev_facets.back(), surface );MB_CHK_SET_ERR( rval, "Failed to resolve the boundary case" ); 01028 } 01029 else 01030 { 01031 // look up nearest facet 01032 01033 // Get OBB Tree for surface 01034 EntityHandle root; 01035 rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's OBB tree root" ); 01036 01037 // Get closest triangle on surface 01038 const CartVect point( xyz ); 01039 CartVect nearest; 01040 EntityHandle facet_out; 01041 rval = geomTopoTool->obb_tree()->closest_to_location( point.array(), root, nearest.array(), facet_out );MB_CHK_SET_ERR( rval, "Failed to find the closest point to location" ); 01042 01043 rval = boundary_case( volume, dir, uvw[0], uvw[1], uvw[2], facet_out, surface );MB_CHK_SET_ERR( rval, "Failed to resolve the boundary case" ); 01044 } 01045 01046 result = dir; 01047 01048 return MB_SUCCESS; 01049 } 01050 01051 // use spherical area test to determine inside/outside of a polyhedron. 01052 ErrorCode GeomQueryTool::point_in_volume_slow( EntityHandle volume, const double xyz[3], int& result ) 01053 { 01054 ErrorCode rval; 01055 Range faces; 01056 std::vector< EntityHandle > surfs; 01057 std::vector< int > senses; 01058 double sum = 0.0; 01059 const CartVect point( xyz ); 01060 01061 rval = MBI->get_child_meshsets( volume, surfs );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" ); 01062 01063 senses.resize( surfs.size() ); 01064 rval = geomTopoTool->get_surface_senses( volume, surfs.size(), &surfs[0], &senses[0] );MB_CHK_SET_ERR( rval, "Failed to get the volume's surface senses" ); 01065 01066 for( unsigned i = 0; i < surfs.size(); ++i ) 01067 { 01068 if( !senses[i] ) // skip non-manifold surfaces 01069 continue; 01070 01071 double surf_area = 0.0, face_area; 01072 faces.clear(); 01073 rval = MBI->get_entities_by_dimension( surfs[i], 2, faces );MB_CHK_SET_ERR( rval, "Failed to get the surface entities by dimension" ); 01074 01075 for( Range::iterator j = faces.begin(); j != faces.end(); ++j ) 01076 { 01077 rval = poly_solid_angle( *j, point, face_area );MB_CHK_SET_ERR( rval, "Failed to determin the polygon's solid angle" ); 01078 01079 surf_area += face_area; 01080 } 01081 01082 sum += senses[i] * surf_area; 01083 } 01084 01085 result = fabs( sum ) > 2.0 * M_PI; 01086 return MB_SUCCESS; 01087 } 01088 01089 ErrorCode GeomQueryTool::find_volume( const double xyz[3], EntityHandle& volume, const double* dir ) 01090 { 01091 ErrorCode rval; 01092 volume = 0; 01093 01094 EntityHandle global_surf_tree_root = geomTopoTool->get_one_vol_root(); 01095 01096 // fast check - make sure point is in the implicit complement bounding box 01097 int ic_result; 01098 EntityHandle ic_handle; 01099 rval = geomTopoTool->get_implicit_complement( ic_handle );MB_CHK_SET_ERR( rval, "Failed to get the implicit complement handle" ); 01100 01101 rval = point_in_box( ic_handle, xyz, ic_result );MB_CHK_SET_ERR( rval, "Failed to check implicit complement for containment" ); 01102 if( ic_result == 0 ) 01103 { 01104 volume = 0; 01105 return MB_ENTITY_NOT_FOUND; 01106 } 01107 01108 // if geomTopoTool doesn't have a global tree, use a loop over vols (slow) 01109 if( !global_surf_tree_root ) 01110 { 01111 rval = find_volume_slow( xyz, volume, dir ); 01112 return rval; 01113 } 01114 01115 moab::CartVect uvw( 0.0 ); 01116 01117 if( dir ) 01118 { 01119 uvw[0] = dir[0]; 01120 uvw[1] = dir[1]; 01121 uvw[2] = dir[2]; 01122 } 01123 01124 if( uvw == 0.0 ) 01125 { 01126 uvw[0] = rand(); 01127 uvw[1] = rand(); 01128 uvw[2] = rand(); 01129 } 01130 01131 // always normalize direction 01132 uvw.normalize(); 01133 01134 // fire a ray along dir and get surface 01135 const double huge_val = std::numeric_limits< double >::max(); 01136 double pos_ray_len = huge_val; 01137 double neg_ray_len = -huge_val; 01138 01139 // RIS output data 01140 std::vector< double > dists; 01141 std::vector< EntityHandle > surfs; 01142 std::vector< EntityHandle > facets; 01143 01144 FindVolumeIntRegCtxt find_vol_reg_ctxt; 01145 OrientedBoxTreeTool::IntersectSearchWindow search_win( &pos_ray_len, &neg_ray_len ); 01146 rval = 01147 geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, global_surf_tree_root, numericalPrecision, 01148 xyz, uvw.array(), search_win, find_vol_reg_ctxt );MB_CHK_SET_ERR( rval, "Failed in global tree ray fire" ); 01149 01150 // if there was no intersection, no volume is found 01151 if( surfs.size() == 0 || surfs[0] == 0 ) 01152 { 01153 volume = 0; 01154 return MB_ENTITY_NOT_FOUND; 01155 } 01156 01157 // get the positive distance facet, surface hit 01158 EntityHandle facet = facets[0]; 01159 EntityHandle surf = surfs[0]; 01160 01161 // get these now, we're going to use them no matter what 01162 EntityHandle fwd_vol, bwd_vol; 01163 rval = geomTopoTool->get_surface_senses( surf, fwd_vol, bwd_vol );MB_CHK_SET_ERR( rval, "Failed to get sense data" ); 01164 EntityHandle parent_vols[2]; 01165 parent_vols[0] = fwd_vol; 01166 parent_vols[1] = bwd_vol; 01167 01168 // get triangle normal 01169 std::vector< EntityHandle > conn; 01170 CartVect coords[3]; 01171 rval = MBI->get_connectivity( &facet, 1, conn );MB_CHK_SET_ERR( rval, "Failed to get triangle connectivity" ); 01172 01173 rval = MBI->get_coords( &conn[0], 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get triangle coordinates" ); 01174 01175 CartVect normal = ( coords[1] - coords[0] ) * ( coords[2] - coords[0] ); 01176 normal.normalize(); 01177 01178 // reverse direction if a hit in the negative direction is found 01179 if( dists[0] < 0 ) { uvw *= -1; } 01180 01181 // if this is a "forward" intersection return the first sense entity 01182 // otherwise return the second, "reverse" sense entity 01183 double dot_prod = uvw % normal; 01184 int idx = dot_prod > 0.0 ? 0 : 1; 01185 01186 if( dot_prod == 0.0 ) 01187 { 01188 std::cerr << "Tangent dot product in find_volume. Shouldn't be here." << std::endl; 01189 volume = 0; 01190 return MB_FAILURE; 01191 } 01192 01193 volume = parent_vols[idx]; 01194 01195 return MB_SUCCESS; 01196 } 01197 01198 ErrorCode GeomQueryTool::find_volume_slow( const double xyz[3], EntityHandle& volume, const double* dir ) 01199 { 01200 ErrorCode rval; 01201 volume = 0; 01202 // get all volumes 01203 Range all_vols; 01204 rval = geomTopoTool->get_gsets_by_dimension( 3, all_vols );MB_CHK_SET_ERR( rval, "Failed to get all volumes in the model" ); 01205 01206 Range::iterator it; 01207 int result = 0; 01208 for( it = all_vols.begin(); it != all_vols.end(); it++ ) 01209 { 01210 rval = point_in_volume( *it, xyz, result, dir );MB_CHK_SET_ERR( rval, "Failed in point in volume loop" ); 01211 if( result ) 01212 { 01213 volume = *it; 01214 break; 01215 } 01216 } 01217 return volume ? MB_SUCCESS : MB_ENTITY_NOT_FOUND; 01218 } 01219 01220 // detemine distance to nearest surface 01221 ErrorCode GeomQueryTool::closest_to_location( EntityHandle volume, const double coords[3], double& result, 01222 EntityHandle* closest_surface ) 01223 { 01224 // Get OBB Tree for volume 01225 EntityHandle root; 01226 ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's obb tree root" ); 01227 01228 // Get closest triangles in volume 01229 const CartVect point( coords ); 01230 CartVect nearest; 01231 EntityHandle facet_out; 01232 01233 rval = geomTopoTool->obb_tree()->closest_to_location( point.array(), root, nearest.array(), facet_out, 01234 closest_surface );MB_CHK_SET_ERR( rval, "Failed to get the closest intersection to location" ); 01235 // calculate distance between point and nearest facet 01236 result = ( point - nearest ).length(); 01237 01238 return MB_SUCCESS; 01239 } 01240 01241 // calculate volume of polyhedron 01242 ErrorCode GeomQueryTool::measure_volume( EntityHandle volume, double& result ) 01243 { 01244 ErrorCode rval; 01245 std::vector< EntityHandle > surfaces; 01246 result = 0.0; 01247 01248 // don't try to calculate volume of implicit complement 01249 if( geomTopoTool->is_implicit_complement( volume ) ) 01250 { 01251 result = 1.0; 01252 return MB_SUCCESS; 01253 } 01254 01255 // get surfaces from volume 01256 rval = MBI->get_child_meshsets( volume, surfaces );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" ); 01257 01258 // get surface senses 01259 std::vector< int > senses( surfaces.size() ); 01260 rval = geomTopoTool->get_surface_senses( volume, surfaces.size(), &surfaces[0], &senses[0] );MB_CHK_SET_ERR( rval, "Failed to retrieve surface-volume sense data. Cannot calculate volume" ); 01261 01262 for( unsigned i = 0; i < surfaces.size(); ++i ) 01263 { 01264 // skip non-manifold surfaces 01265 if( !senses[i] ) continue; 01266 01267 // get triangles in surface 01268 Range triangles; 01269 rval = MBI->get_entities_by_dimension( surfaces[i], 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" ); 01270 01271 if( !triangles.all_of_type( MBTRI ) ) 01272 { 01273 std::cout << "WARNING: Surface " << surfaces[i] // todo: use geomtopotool to get id by entity handle 01274 << " contains non-triangle elements. Volume calculation may be incorrect." << std::endl; 01275 triangles.clear(); 01276 rval = MBI->get_entities_by_type( surfaces[i], MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" ); 01277 } 01278 01279 // calculate signed volume beneath surface (x 6.0) 01280 double surf_sum = 0.0; 01281 const EntityHandle* conn; 01282 int len; 01283 CartVect coords[3]; 01284 for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j ) 01285 { 01286 rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the current triangle" ); 01287 if( 3 != len ) { MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" ); } 01288 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get the coordinates of the current triangle's vertices" ); 01289 01290 coords[1] -= coords[0]; 01291 coords[2] -= coords[0]; 01292 surf_sum += ( coords[0] % ( coords[1] * coords[2] ) ); 01293 } 01294 result += senses[i] * surf_sum; 01295 } 01296 01297 result /= 6.0; 01298 return MB_SUCCESS; 01299 } 01300 01301 // sum area of elements in surface 01302 ErrorCode GeomQueryTool::measure_area( EntityHandle surface, double& result ) 01303 { 01304 // get triangles in surface 01305 Range triangles; 01306 ErrorCode rval = MBI->get_entities_by_dimension( surface, 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface entities" ); 01307 if( !triangles.all_of_type( MBTRI ) ) 01308 { 01309 std::cout << "WARNING: Surface " << surface // todo: use geomtopotool to get id by entity handle 01310 << " contains non-triangle elements. Area calculation may be incorrect." << std::endl; 01311 triangles.clear(); 01312 rval = MBI->get_entities_by_type( surface, MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to the surface's triangle entities" ); 01313 } 01314 01315 // calculate sum of area of triangles 01316 result = 0.0; 01317 const EntityHandle* conn; 01318 int len; 01319 CartVect coords[3]; 01320 for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j ) 01321 { 01322 rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's connectivity" ); 01323 if( 3 != len ) { MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" ); } 01324 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's vertex coordinates" ); 01325 01326 // calculated area using cross product of triangle edges 01327 CartVect v1 = coords[1] - coords[0]; 01328 CartVect v2 = coords[2] - coords[0]; 01329 CartVect xp = v1 * v2; 01330 result += xp.length(); 01331 } 01332 result *= 0.5; 01333 return MB_SUCCESS; 01334 } 01335 01336 ErrorCode GeomQueryTool::get_normal( EntityHandle surf, const double in_pt[3], double angle[3], 01337 const RayHistory* history ) 01338 { 01339 EntityHandle root; 01340 ErrorCode rval = geomTopoTool->get_root( surf, root );MB_CHK_SET_ERR( rval, "Failed to get the surface's obb tree root" ); 01341 01342 std::vector< EntityHandle > facets; 01343 01344 // if no history or history empty, use nearby facets 01345 if( !history || ( history->prev_facets.size() == 0 ) ) 01346 { 01347 rval = geomTopoTool->obb_tree()->closest_to_location( in_pt, root, numericalPrecision, facets );MB_CHK_SET_ERR( rval, "Failed to get closest intersection to location" ); 01348 } 01349 // otherwise use most recent facet in history 01350 else 01351 { 01352 facets.push_back( history->prev_facets.back() ); 01353 } 01354 01355 CartVect coords[3], normal( 0.0 ); 01356 const EntityHandle* conn; 01357 int len; 01358 for( unsigned i = 0; i < facets.size(); ++i ) 01359 { 01360 rval = MBI->get_connectivity( facets[i], conn, len );MB_CHK_SET_ERR( rval, "Failed to get facet connectivity" ); 01361 if( 3 != len ) { MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" ); } 01362 01363 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" ); 01364 01365 coords[1] -= coords[0]; 01366 coords[2] -= coords[0]; 01367 normal += coords[1] * coords[2]; 01368 } 01369 01370 normal.normalize(); 01371 normal.get( angle ); 01372 01373 return MB_SUCCESS; 01374 } 01375 01376 /* SECTION II (private) */ 01377 01378 // If point is on boundary, then this function is called to 01379 // discriminate cases in which the ray is entering or leaving. 01380 // result= 1 -> inside volume or entering volume 01381 // result= 0 -> outside volume or leaving volume 01382 // result=-1 -> on boundary with null or tangent uvw 01383 ErrorCode GeomQueryTool::boundary_case( EntityHandle volume, int& result, double u, double v, double w, 01384 EntityHandle facet, EntityHandle surface ) 01385 { 01386 ErrorCode rval; 01387 01388 // test to see if uvw is provided 01389 if( u <= 1.0 && v <= 1.0 && w <= 1.0 ) 01390 { 01391 01392 const CartVect ray_vector( u, v, w ); 01393 CartVect coords[3], normal( 0.0 ); 01394 const EntityHandle* conn; 01395 int len, sense_out; 01396 01397 rval = MBI->get_connectivity( facet, conn, len );MB_CHK_SET_ERR( rval, "Failed to get the triangle's connectivity" ); 01398 if( 3 != len ) { MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" ); } 01399 01400 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" ); 01401 01402 rval = geomTopoTool->get_sense( surface, volume, sense_out );MB_CHK_SET_ERR( rval, "Failed to get the surface's sense with respect to it's volume" ); 01403 01404 coords[1] -= coords[0]; 01405 coords[2] -= coords[0]; 01406 normal = sense_out * ( coords[1] * coords[2] ); 01407 01408 double sense = ray_vector % normal; 01409 01410 if( sense < 0.0 ) 01411 { 01412 result = 1; // inside or entering 01413 } 01414 else if( sense > 0.0 ) 01415 { 01416 result = 0; // outside or leaving 01417 } 01418 else if( sense == 0.0 ) 01419 { 01420 result = -1; // tangent, therefore on boundary 01421 } 01422 else 01423 { 01424 result = -1; // failure 01425 MB_SET_ERR( MB_FAILURE, "Failed to resolve boundary case" ); 01426 } 01427 01428 // if uvw not provided, return on_boundary. 01429 } 01430 else 01431 { 01432 result = -1; // on boundary 01433 return MB_SUCCESS; 01434 } 01435 01436 return MB_SUCCESS; 01437 } 01438 01439 // point_in_volume_slow, including poly_solid_angle helper subroutine 01440 // are adapted from "Point in Polyhedron Testing Using Spherical Polygons", Paulo Cezar 01441 // Pinto Carvalho and Paulo Roma Cavalcanti, _Graphics Gems V_, pg. 42. Original algorithm 01442 // was described in "An Efficient Point In Polyhedron Algorithm", Jeff Lane, Bob Magedson, 01443 // and Mike Rarick, _Computer Vision, Graphics, and Image Processing 26_, pg. 118-225, 1984. 01444 01445 // helper function for point_in_volume_slow. calculate area of a polygon 01446 // projected into a unit-sphere space 01447 ErrorCode GeomQueryTool::poly_solid_angle( EntityHandle face, const CartVect& point, double& area ) 01448 { 01449 ErrorCode rval; 01450 01451 // Get connectivity 01452 const EntityHandle* conn; 01453 int len; 01454 rval = MBI->get_connectivity( face, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the polygon" ); 01455 01456 // Allocate space to store vertices 01457 CartVect coords_static[4]; 01458 std::vector< CartVect > coords_dynamic; 01459 CartVect* coords = coords_static; 01460 if( (unsigned)len > ( sizeof( coords_static ) / sizeof( coords_static[0] ) ) ) 01461 { 01462 coords_dynamic.resize( len ); 01463 coords = &coords_dynamic[0]; 01464 } 01465 01466 // get coordinates 01467 rval = MBI->get_coords( conn, len, coords->array() );MB_CHK_SET_ERR( rval, "Failed to get the coordinates of the polygon vertices" ); 01468 01469 // calculate normal 01470 CartVect norm( 0.0 ), v1, v0 = coords[1] - coords[0]; 01471 for( int i = 2; i < len; ++i ) 01472 { 01473 v1 = coords[i] - coords[0]; 01474 norm += v0 * v1; 01475 v0 = v1; 01476 } 01477 01478 // calculate area 01479 double s, ang; 01480 area = 0.0; 01481 CartVect r, n1, n2, b, a = coords[len - 1] - coords[0]; 01482 for( int i = 0; i < len; ++i ) 01483 { 01484 r = coords[i] - point; 01485 b = coords[( i + 1 ) % len] - coords[i]; 01486 n1 = a * r; // = norm1 (magnitude is important) 01487 n2 = r * b; // = norm2 (magnitude is important) 01488 s = ( n1 % n2 ) / ( n1.length() * n2.length() ); // = cos(angle between norm1,norm2) 01489 ang = s <= -1.0 ? M_PI : s >= 1.0 ? 0.0 : acos( s ); // = acos(s) 01490 s = ( b * a ) % norm; // =orientation of triangle wrt point 01491 area += s > 0.0 ? M_PI - ang : M_PI + ang; 01492 a = -b; 01493 } 01494 01495 area -= M_PI * ( len - 2 ); 01496 if( ( norm % r ) > 0 ) area = -area; 01497 return MB_SUCCESS; 01498 } 01499 01500 void GeomQueryTool::set_overlap_thickness( double new_thickness ) 01501 { 01502 01503 if( new_thickness < 0 || new_thickness > 100 ) 01504 { std::cerr << "Invalid overlap_thickness = " << new_thickness << std::endl; } 01505 else 01506 { 01507 overlapThickness = new_thickness; 01508 } 01509 std::cout << "Set overlap thickness = " << overlapThickness << std::endl; 01510 } 01511 01512 void GeomQueryTool::set_numerical_precision( double new_precision ) 01513 { 01514 01515 if( new_precision <= 0 || new_precision > 1 ) 01516 { std::cerr << "Invalid numerical_precision = " << numericalPrecision << std::endl; } 01517 else 01518 { 01519 numericalPrecision = new_precision; 01520 } 01521 01522 std::cout << "Set numerical precision = " << numericalPrecision << std::endl; 01523 } 01524 01525 } // namespace moab