 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
Tool for querying different aspects of geometric topology sets in MOAB. More...
#include <GeomQueryTool.hpp>
Collaboration diagram for moab::GeomQueryTool:Classes | |
| class | RayHistory |
Public Member Functions | |
| GeomQueryTool (Interface *impl, bool find_geoments=false, EntityHandle modelRootSet=0, bool p_rootSets_vector=true, bool restore_rootSets=true, bool trace_counting=false, double overlap_thickness=0., double numerical_precision=0.001) | |
| GeomQueryTool (GeomTopoTool *geomtopotool, bool trace_counting=false, double overlap_thickness=0., double numerical_precision=0.001) | |
| ~GeomQueryTool () | |
| ErrorCode | initialize () |
| ErrorCode | ray_fire (const EntityHandle volume, const double ray_start[3], const double ray_dir[3], EntityHandle &next_surf, double &next_surf_dist, RayHistory *history=NULL, double dist_limit=0, int ray_orientation=1, OrientedBoxTreeTool::TrvStats *stats=NULL) |
| find the next surface crossing from a given point in a given direction More... | |
| ErrorCode | point_in_volume (const EntityHandle volume, const double xyz[3], int &result, const double *uvw=NULL, const RayHistory *history=NULL) |
| Test if a point is inside or outside a volume. More... | |
| ErrorCode | point_in_volume_slow (const EntityHandle volume, const double xyz[3], int &result) |
| Robust test if a point is inside or outside a volume using unit sphere area method. More... | |
| ErrorCode | find_volume (const double xyz[3], EntityHandle &volume, const double *dir=NULL) |
| Find volume for a given location. More... | |
| ErrorCode | find_volume_slow (const double xyz[3], EntityHandle &volume, const double *dir=NULL) |
| Find volume for a given location using loop. (slow) More... | |
| ErrorCode | point_in_box (const EntityHandle volume, const double point[3], int &inside) |
| Test if a point is inside or outsize a volume's axis-aligned bounding box. More... | |
| ErrorCode | test_volume_boundary (const EntityHandle volume, const EntityHandle surface, const double xyz[3], const double uvw[3], int &result, const RayHistory *history=NULL) |
| Given a ray starting at a surface of a volume, check whether the ray enters or exits the volume. More... | |
| ErrorCode | closest_to_location (EntityHandle volume, const double point[3], double &result, EntityHandle *closest_surface=0) |
| Find the distance to the point on the boundary of the volume closest to the test point. More... | |
| ErrorCode | measure_volume (EntityHandle volume, double &result) |
| ErrorCode | measure_area (EntityHandle surface, double &result) |
| ErrorCode | get_normal (EntityHandle surf, const double xyz[3], double angle[3], const RayHistory *history=NULL) |
| void | set_overlap_thickness (double new_overlap_thickness) |
| State object used in calls to ray_fire() More... | |
| void | set_numerical_precision (double new_precision) |
| void | set_verbosity (bool value) |
| bool | get_verbosity () const |
| double | get_numerical_precision () |
| double | get_overlap_thickness () |
| GeomTopoTool * | gttool () |
| Interface * | moab_instance () |
Private Member Functions | |
| ErrorCode | boundary_case (EntityHandle volume, int &result, double u, double v, double w, EntityHandle facet, EntityHandle surface) |
| determine the point membership when the point is effectively on the boundary More... | |
| ErrorCode | poly_solid_angle (EntityHandle face, const CartVect &point, double &area) |
Private Attributes | |
| GeomTopoTool * | geomTopoTool |
| bool | verbose |
| bool | owns_gtt |
| Interface * | MBI |
| OrientedBoxTreeTool * | obbTreeTool |
| bool | counting |
| long long int | n_pt_in_vol_calls |
| long long int | n_ray_fire_calls |
| double | overlapThickness |
| double | numericalPrecision |
| Tag | senseTag |
Tool for querying different aspects of geometric topology sets in MOAB.
Given the conventions established in GeomTopoTool for representing geometric topology through a hierarchy of meshsets, this tool provides a set of methods to query different geometric aspects of those geometric topology sets including:
A feature of these queries is that there is some support for overlapping geometries.
Definition at line 45 of file GeomQueryTool.hpp.
| moab::GeomQueryTool::GeomQueryTool | ( | Interface * | impl, |
| bool | find_geoments = false, |
||
| EntityHandle | modelRootSet = 0, |
||
| bool | p_rootSets_vector = true, |
||
| bool | restore_rootSets = true, |
||
| bool | trace_counting = false, |
||
| double | overlap_thickness = 0., |
||
| double | numerical_precision = 0.001 |
||
| ) |
Definition at line 624 of file GeomQueryTool.cpp.
632 : verbose( false ), owns_gtt( true )
633 {
634 geomTopoTool = new GeomTopoTool( impl, find_geomsets, modelRootSet, p_rootSets_vector, restore_rootSets );
635
636 senseTag = geomTopoTool->get_sense_tag();
637
638 obbTreeTool = geomTopoTool->obb_tree();
639 MBI = geomTopoTool->get_moab_instance();
640
641 counting = trace_counting;
642 overlapThickness = overlap_thickness;
643 numericalPrecision = numerical_precision;
644
645 // reset query counters
646 n_pt_in_vol_calls = 0;
647 n_ray_fire_calls = 0;
648 }
References counting, geomTopoTool, moab::GeomTopoTool::get_moab_instance(), moab::GeomTopoTool::get_sense_tag(), MBI, n_pt_in_vol_calls, n_ray_fire_calls, numericalPrecision, moab::GeomTopoTool::obb_tree(), obbTreeTool, overlapThickness, and senseTag.
| moab::GeomQueryTool::GeomQueryTool | ( | GeomTopoTool * | geomtopotool, |
| bool | trace_counting = false, |
||
| double | overlap_thickness = 0., |
||
| double | numerical_precision = 0.001 |
||
| ) |
Definition at line 650 of file GeomQueryTool.cpp.
654 : verbose( false ), owns_gtt( false )
655 {
656
657 geomTopoTool = geomtopotool;
658
659 senseTag = geomTopoTool->get_sense_tag();
660
661 obbTreeTool = geomTopoTool->obb_tree();
662 MBI = geomTopoTool->get_moab_instance();
663
664 counting = trace_counting;
665 overlapThickness = overlap_thickness;
666 numericalPrecision = numerical_precision;
667
668 // reset query counters
669 n_pt_in_vol_calls = 0;
670 n_ray_fire_calls = 0;
671 }
References counting, geomTopoTool, moab::GeomTopoTool::get_moab_instance(), moab::GeomTopoTool::get_sense_tag(), MBI, n_pt_in_vol_calls, n_ray_fire_calls, numericalPrecision, moab::GeomTopoTool::obb_tree(), obbTreeTool, overlapThickness, and senseTag.
| moab::GeomQueryTool::~GeomQueryTool | ( | ) |
Definition at line 673 of file GeomQueryTool.cpp.
674 {
675 if( owns_gtt )
676 {
677 delete geomTopoTool;
678 }
679 }
References geomTopoTool, and owns_gtt.
|
private |
determine the point membership when the point is effectively on the boundary
Called by point_in_volume when the point is with tolerance of the boundary. Compares the ray direction with the surface normal to determine a volume membership.
Definition at line 1488 of file GeomQueryTool.cpp.
1495 {
1496 ErrorCode rval;
1497
1498 // test to see if uvw is provided
1499 if( u <= 1.0 && v <= 1.0 && w <= 1.0 )
1500 {
1501
1502 const CartVect ray_vector( u, v, w );
1503 CartVect coords[3], normal( 0.0 );
1504 const EntityHandle* conn;
1505 int len, sense_out;
1506
1507 rval = MBI->get_connectivity( facet, conn, len );MB_CHK_SET_ERR( rval, "Failed to get the triangle's connectivity" );
1508 if( 3 != len )
1509 {
1510 MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
1511 }
1512
1513 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" );
1514
1515 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" );
1516
1517 coords[1] -= coords[0];
1518 coords[2] -= coords[0];
1519 normal = sense_out * ( coords[1] * coords[2] );
1520
1521 double sense = ray_vector % normal;
1522
1523 if( sense < 0.0 )
1524 {
1525 result = 1; // inside or entering
1526 }
1527 else if( sense > 0.0 )
1528 {
1529 result = 0; // outside or leaving
1530 }
1531 else if( sense == 0.0 )
1532 {
1533 result = -1; // tangent, therefore on boundary
1534 }
1535 else
1536 {
1537 result = -1; // failure
1538 MB_SET_ERR( MB_FAILURE, "Failed to resolve boundary case" );
1539 }
1540
1541 // if uvw not provided, return on_boundary.
1542 }
1543 else
1544 {
1545 result = -1; // on boundary
1546 return MB_SUCCESS;
1547 }
1548
1549 return MB_SUCCESS;
1550 }
References ErrorCode, geomTopoTool, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::GeomTopoTool::get_sense(), MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, and MBI.
Referenced by point_in_volume(), and test_volume_boundary().
| ErrorCode moab::GeomQueryTool::closest_to_location | ( | EntityHandle | volume, |
| const double | point[3], | ||
| double & | result, | ||
| EntityHandle * | closest_surface = 0 |
||
| ) |
Find the distance to the point on the boundary of the volume closest to the test point.
| volume | Volume to query |
| point | Coordinates of test point |
| result | Set to the minimum distance from point to a surface in volume |
Definition at line 1313 of file GeomQueryTool.cpp.
1317 {
1318 // Get OBB Tree for volume
1319 EntityHandle root;
1320 ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's obb tree root" );
1321
1322 // Get closest triangles in volume
1323 const CartVect point( coords );
1324 CartVect nearest;
1325 EntityHandle facet_out;
1326
1327 rval = geomTopoTool->obb_tree()->closest_to_location( point.array(), root, nearest.array(), facet_out,
1328 closest_surface );MB_CHK_SET_ERR( rval, "Failed to get the closest intersection to location" );
1329 // calculate distance between point and nearest facet
1330 result = ( point - nearest ).length();
1331
1332 return MB_SUCCESS;
1333 }
References moab::CartVect::array(), moab::OrientedBoxTreeTool::closest_to_location(), ErrorCode, geomTopoTool, moab::GeomTopoTool::get_root(), length(), MB_CHK_SET_ERR, MB_SUCCESS, and moab::GeomTopoTool::obb_tree().
| ErrorCode moab::GeomQueryTool::find_volume | ( | const double | xyz[3], |
| EntityHandle & | volume, | ||
| const double * | dir = NULL |
||
| ) |
Find volume for a given location.
Determines which volume contains the point if possible. If no volume is found, a null EntityHandle is returned along with a MB_ENTITY_NOT_FOUND ErrorCode.
| xyz | The location to test |
| volume | Set to volume handle containing the location if found |
| dir | Optional direction to use for underlying ray fire query. Used to ensure consistent results when a ray direction is known. If NULL or {0,0,0} is given, a random direction will be used. |
Definition at line 1178 of file GeomQueryTool.cpp.
1179 {
1180 ErrorCode rval;
1181 volume = 0;
1182
1183 EntityHandle global_surf_tree_root = geomTopoTool->get_one_vol_root();
1184
1185 // fast check - make sure point is in the implicit complement bounding box
1186 int ic_result;
1187 EntityHandle ic_handle;
1188 rval = geomTopoTool->get_implicit_complement( ic_handle );MB_CHK_SET_ERR( rval, "Failed to get the implicit complement handle" );
1189
1190 rval = point_in_box( ic_handle, xyz, ic_result );MB_CHK_SET_ERR( rval, "Failed to check implicit complement for containment" );
1191 if( ic_result == 0 )
1192 {
1193 volume = 0;
1194 return MB_ENTITY_NOT_FOUND;
1195 }
1196
1197 // if geomTopoTool doesn't have a global tree, use a loop over vols (slow)
1198 if( !global_surf_tree_root )
1199 {
1200 rval = find_volume_slow( xyz, volume, dir );
1201 return rval;
1202 }
1203
1204 moab::CartVect uvw( 0.0 );
1205
1206 if( dir )
1207 {
1208 uvw[0] = dir[0];
1209 uvw[1] = dir[1];
1210 uvw[2] = dir[2];
1211 }
1212
1213 if( uvw == 0.0 )
1214 {
1215 uvw[0] = rand();
1216 uvw[1] = rand();
1217 uvw[2] = rand();
1218 }
1219
1220 // always normalize direction
1221 uvw.normalize();
1222
1223 // fire a ray along dir and get surface
1224 const double huge_val = std::numeric_limits< double >::max();
1225 double pos_ray_len = huge_val;
1226 double neg_ray_len = -huge_val;
1227
1228 // RIS output data
1229 std::vector< double > dists;
1230 std::vector< EntityHandle > surfs;
1231 std::vector< EntityHandle > facets;
1232
1233 FindVolumeIntRegCtxt find_vol_reg_ctxt;
1234 OrientedBoxTreeTool::IntersectSearchWindow search_win( &pos_ray_len, &neg_ray_len );
1235 rval =
1236 geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, global_surf_tree_root, numericalPrecision,
1237 xyz, uvw.array(), search_win, find_vol_reg_ctxt );MB_CHK_SET_ERR( rval, "Failed in global tree ray fire" );
1238
1239 // if there was no intersection, no volume is found
1240 if( surfs.size() == 0 || surfs[0] == 0 )
1241 {
1242 volume = 0;
1243 return MB_ENTITY_NOT_FOUND;
1244 }
1245
1246 // get the positive distance facet, surface hit
1247 EntityHandle facet = facets[0];
1248 EntityHandle surf = surfs[0];
1249
1250 // get these now, we're going to use them no matter what
1251 EntityHandle fwd_vol, bwd_vol;
1252 rval = geomTopoTool->get_surface_senses( surf, fwd_vol, bwd_vol );MB_CHK_SET_ERR( rval, "Failed to get sense data" );
1253 EntityHandle parent_vols[2];
1254 parent_vols[0] = fwd_vol;
1255 parent_vols[1] = bwd_vol;
1256
1257 // get triangle normal
1258 std::vector< EntityHandle > conn;
1259 CartVect coords[3];
1260 rval = MBI->get_connectivity( &facet, 1, conn );MB_CHK_SET_ERR( rval, "Failed to get triangle connectivity" );
1261
1262 rval = MBI->get_coords( &conn[0], 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get triangle coordinates" );
1263
1264 CartVect normal = ( coords[1] - coords[0] ) * ( coords[2] - coords[0] );
1265 normal.normalize();
1266
1267 // reverse direction if a hit in the negative direction is found
1268 if( dists[0] < 0 )
1269 {
1270 uvw *= -1;
1271 }
1272
1273 // if this is a "forward" intersection return the first sense entity
1274 // otherwise return the second, "reverse" sense entity
1275 double dot_prod = uvw % normal;
1276 int idx = dot_prod > 0.0 ? 0 : 1;
1277
1278 if( dot_prod == 0.0 )
1279 {
1280 std::cerr << "Tangent dot product in find_volume. Shouldn't be here." << std::endl;
1281 volume = 0;
1282 return MB_FAILURE;
1283 }
1284
1285 volume = parent_vols[idx];
1286
1287 return MB_SUCCESS;
1288 }
References moab::CartVect::array(), ErrorCode, find_volume_slow(), geomTopoTool, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::GeomTopoTool::get_implicit_complement(), moab::GeomTopoTool::get_one_vol_root(), moab::GeomTopoTool::get_surface_senses(), MB_CHK_SET_ERR, MB_ENTITY_NOT_FOUND, MB_SUCCESS, MBI, moab::CartVect::normalize(), numericalPrecision, moab::GeomTopoTool::obb_tree(), point_in_box(), and moab::OrientedBoxTreeTool::ray_intersect_sets().
| ErrorCode moab::GeomQueryTool::find_volume_slow | ( | const double | xyz[3], |
| EntityHandle & | volume, | ||
| const double * | dir = NULL |
||
| ) |
Find volume for a given location using loop. (slow)
Loops over all volumes in the model, checking for point containment
| xyz | The location to test |
| volume | Set to volume handle containing the location if found |
| dir | Optional direction to use for underlying ray fire query. Used to ensure consistent results when a ray direction is known. If NULL or {0,0,0} is given, a random direction will be used. |
Definition at line 1290 of file GeomQueryTool.cpp.
1291 {
1292 ErrorCode rval;
1293 volume = 0;
1294 // get all volumes
1295 Range all_vols;
1296 rval = geomTopoTool->get_gsets_by_dimension( 3, all_vols );MB_CHK_SET_ERR( rval, "Failed to get all volumes in the model" );
1297
1298 Range::iterator it;
1299 int result = 0;
1300 for( it = all_vols.begin(); it != all_vols.end(); it++ )
1301 {
1302 rval = point_in_volume( *it, xyz, result, dir );MB_CHK_SET_ERR( rval, "Failed in point in volume loop" );
1303 if( result )
1304 {
1305 volume = *it;
1306 break;
1307 }
1308 }
1309 return volume ? MB_SUCCESS : MB_ENTITY_NOT_FOUND;
1310 }
References moab::Range::begin(), moab::Range::end(), ErrorCode, geomTopoTool, moab::GeomTopoTool::get_gsets_by_dimension(), MB_CHK_SET_ERR, MB_ENTITY_NOT_FOUND, MB_SUCCESS, and point_in_volume().
Referenced by find_volume().
| ErrorCode moab::GeomQueryTool::get_normal | ( | EntityHandle | surf, |
| const double | xyz[3], | ||
| double | angle[3], | ||
| const RayHistory * | history = NULL |
||
| ) |
Get the normal to a given surface at the point on the surface closest to a given point
This method first identifies which facet contains this point and then calculates the unit outward normal of that facet. The facet of the provided volume that is nearest the provided point is used for this calculation. The search for that facet can be circumvented by providing a RayHistory, in which case the last facet of the history will be used.
| surf | Surface on which to get normal |
| xyz | Point on surf |
| angle | Set to coordinates of surface normal nearest xyz |
| history | Optional ray history from a previous call to ray_fire(). If present and non-empty, return the normal of the most recently intersected facet, ignoring xyz. |
Definition at line 1436 of file GeomQueryTool.cpp.
1440 {
1441 EntityHandle root;
1442 ErrorCode rval = geomTopoTool->get_root( surf, root );MB_CHK_SET_ERR( rval, "Failed to get the surface's obb tree root" );
1443
1444 std::vector< EntityHandle > facets;
1445
1446 // if no history or history empty, use nearby facets
1447 if( !history || ( history->prev_facets.size() == 0 ) )
1448 {
1449 rval = geomTopoTool->obb_tree()->closest_to_location( in_pt, root, numericalPrecision, facets );MB_CHK_SET_ERR( rval, "Failed to get closest intersection to location" );
1450 }
1451 // otherwise use most recent facet in history
1452 else
1453 {
1454 facets.push_back( history->prev_facets.back() );
1455 }
1456
1457 CartVect coords[3], normal( 0.0 );
1458 const EntityHandle* conn;
1459 int len;
1460 for( unsigned i = 0; i < facets.size(); ++i )
1461 {
1462 rval = MBI->get_connectivity( facets[i], conn, len );MB_CHK_SET_ERR( rval, "Failed to get facet connectivity" );
1463 if( 3 != len )
1464 {
1465 MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
1466 }
1467
1468 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" );
1469
1470 coords[1] -= coords[0];
1471 coords[2] -= coords[0];
1472 normal += coords[1] * coords[2];
1473 }
1474
1475 normal.normalize();
1476 normal.get( angle );
1477
1478 return MB_SUCCESS;
1479 }
References moab::angle(), moab::OrientedBoxTreeTool::closest_to_location(), ErrorCode, geomTopoTool, moab::CartVect::get(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::GeomTopoTool::get_root(), MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, MBI, moab::CartVect::normalize(), numericalPrecision, moab::GeomTopoTool::obb_tree(), and moab::GeomQueryTool::RayHistory::prev_facets.
|
inline |
Definition at line 376 of file GeomQueryTool.hpp.
377 {
378 return numericalPrecision;
379 }
References numericalPrecision.
|
inline |
Definition at line 381 of file GeomQueryTool.hpp.
382 {
383 return overlapThickness;
384 }
References overlapThickness.
|
inline |
Definition at line 371 of file GeomQueryTool.hpp.
372 {
373 return verbose;
374 }
References verbose.
|
inline |
Definition at line 386 of file GeomQueryTool.hpp.
387 {
388 return geomTopoTool;
389 }
References geomTopoTool.
| ErrorCode moab::GeomQueryTool::initialize | ( | ) |
Definition at line 681 of file GeomQueryTool.cpp.
682 {
683
684 ErrorCode rval;
685
686 rval = geomTopoTool->find_geomsets();MB_CHK_SET_ERR( rval, "Failed to find geometry sets" );
687
688 rval = geomTopoTool->setup_implicit_complement();MB_CHK_SET_ERR( rval, "Couldn't setup the implicit complement" );
689
690 rval = geomTopoTool->construct_obb_trees();MB_CHK_SET_ERR( rval, "Failed to construct OBB trees" );
691
692 return MB_SUCCESS;
693 }
References moab::GeomTopoTool::construct_obb_trees(), ErrorCode, moab::GeomTopoTool::find_geomsets(), geomTopoTool, MB_CHK_SET_ERR, MB_SUCCESS, and moab::GeomTopoTool::setup_implicit_complement().
| ErrorCode moab::GeomQueryTool::measure_area | ( | EntityHandle | surface, |
| double & | result | ||
| ) |
Calculate sum of area of triangles
Definition at line 1399 of file GeomQueryTool.cpp.
1400 {
1401 // get triangles in surface
1402 Range triangles;
1403 ErrorCode rval = MBI->get_entities_by_dimension( surface, 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface entities" );
1404 if( !triangles.all_of_type( MBTRI ) )
1405 {
1406 std::cout << "WARNING: Surface " << surface // todo: use geomtopotool to get id by entity handle
1407 << " contains non-triangle elements. Area calculation may be incorrect." << std::endl;
1408 triangles.clear();
1409 rval = MBI->get_entities_by_type( surface, MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to the surface's triangle entities" );
1410 }
1411
1412 // calculate sum of area of triangles
1413 result = 0.0;
1414 const EntityHandle* conn;
1415 int len;
1416 CartVect coords[3];
1417 for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j )
1418 {
1419 rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's connectivity" );
1420 if( 3 != len )
1421 {
1422 MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
1423 }
1424 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's vertex coordinates" );
1425
1426 // calculated area using cross product of triangle edges
1427 CartVect v1 = coords[1] - coords[0];
1428 CartVect v2 = coords[2] - coords[0];
1429 CartVect xp = v1 * v2;
1430 result += xp.length();
1431 }
1432 result *= 0.5;
1433 return MB_SUCCESS;
1434 }
References moab::Range::all_of_type(), moab::Range::begin(), moab::Range::clear(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::get_entities_by_type(), moab::CartVect::length(), MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, MBI, and MBTRI.
| ErrorCode moab::GeomQueryTool::measure_volume | ( | EntityHandle | volume, |
| double & | result | ||
| ) |
Calculate the volume contained in a 'volume'
Definition at line 1336 of file GeomQueryTool.cpp.
1337 {
1338 ErrorCode rval;
1339 std::vector< EntityHandle > surfaces;
1340 result = 0.0;
1341
1342 // don't try to calculate volume of implicit complement
1343 if( geomTopoTool->is_implicit_complement( volume ) )
1344 {
1345 result = 1.0;
1346 return MB_SUCCESS;
1347 }
1348
1349 // get surfaces from volume
1350 rval = MBI->get_child_meshsets( volume, surfaces );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" );
1351
1352 // get surface senses
1353 std::vector< int > senses( surfaces.size() );
1354 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" );
1355
1356 for( unsigned i = 0; i < surfaces.size(); ++i )
1357 {
1358 // skip non-manifold surfaces
1359 if( !senses[i] ) continue;
1360
1361 // get triangles in surface
1362 Range triangles;
1363 rval = MBI->get_entities_by_dimension( surfaces[i], 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" );
1364
1365 if( !triangles.all_of_type( MBTRI ) )
1366 {
1367 std::cout << "WARNING: Surface " << surfaces[i] // todo: use geomtopotool to get id by entity handle
1368 << " contains non-triangle elements. Volume calculation may be incorrect." << std::endl;
1369 triangles.clear();
1370 rval = MBI->get_entities_by_type( surfaces[i], MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" );
1371 }
1372
1373 // calculate signed volume beneath surface (x 6.0)
1374 double surf_sum = 0.0;
1375 const EntityHandle* conn;
1376 int len;
1377 CartVect coords[3];
1378 for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j )
1379 {
1380 rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the current triangle" );
1381 if( 3 != len )
1382 {
1383 MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
1384 }
1385 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" );
1386
1387 coords[1] -= coords[0];
1388 coords[2] -= coords[0];
1389 surf_sum += ( coords[0] % ( coords[1] * coords[2] ) );
1390 }
1391 result += senses[i] * surf_sum;
1392 }
1393
1394 result /= 6.0;
1395 return MB_SUCCESS;
1396 }
References moab::Range::all_of_type(), moab::Range::begin(), moab::Range::clear(), moab::Range::end(), ErrorCode, geomTopoTool, moab::Interface::get_child_meshsets(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::get_entities_by_type(), moab::GeomTopoTool::get_surface_senses(), moab::GeomTopoTool::is_implicit_complement(), MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, MBI, and MBTRI.
|
inline |
| ErrorCode moab::GeomQueryTool::point_in_box | ( | const EntityHandle | volume, |
| const double | point[3], | ||
| int & | inside | ||
| ) |
Test if a point is inside or outsize a volume's axis-aligned bounding box.
For the volume pointed to and the point wished to be tested, returns whether the point is inside or outside the bounding box of the volume. inside = 0, not inside, inside = 1, inside.
This is used as a fast way to determine whether a point is inside or outside of a volume before performing a point_in_volume test which involves firing a ray.
| volume | The volume to test |
| point | The location to test for bounding box containment |
| inside | set to 0 if point is outside the box, 1 if inside |
Definition at line 1077 of file GeomQueryTool.cpp.
1078 {
1079 double minpt[3];
1080 double maxpt[3];
1081 ErrorCode rval = geomTopoTool->get_bounding_coords( volume, minpt, maxpt );MB_CHK_SET_ERR( rval, "Failed to get the bounding coordinates of the volume" );
1082
1083 // early exits
1084 if( point[0] > maxpt[0] || point[0] < minpt[0] )
1085 {
1086 inside = 0;
1087 return rval;
1088 }
1089 if( point[1] > maxpt[1] || point[1] < minpt[1] )
1090 {
1091 inside = 0;
1092 return rval;
1093 }
1094 if( point[2] > maxpt[2] || point[2] < minpt[2] )
1095 {
1096 inside = 0;
1097 return rval;
1098 }
1099 inside = 1;
1100 return rval;
1101 }
References ErrorCode, geomTopoTool, moab::GeomTopoTool::get_bounding_coords(), and MB_CHK_SET_ERR.
Referenced by find_volume(), and point_in_volume().
| ErrorCode moab::GeomQueryTool::point_in_volume | ( | const EntityHandle | volume, |
| const double | xyz[3], | ||
| int & | result, | ||
| const double * | uvw = NULL, |
||
| const RayHistory * | history = NULL |
||
| ) |
Test if a point is inside or outside a volume.
This method finds the point on the boundary of the volume that is nearest the test point (x,y,z). If that point is "close" to a surface, a boundary test is performed based on the normal of the surface at that point and the optional ray direction (u,v,w).
| volume | The volume to test |
| xyz | The location to test for volume containment |
| result | Set to 0 if xyz it outside volume, 1 if inside, and -1 if on boundary. |
| Optional | direction to use for underlying ray fire query. Used to ensure consistent results when a ray direction is known. If NULL or {0,0,0} is given, a random direction will be used. |
| history | Optional RayHistory object to pass to underlying ray fire query. The history is not modified by this call. |
Definition at line 915 of file GeomQueryTool.cpp.
920 {
921 // take some stats that are independent of nps
922 if( counting ) ++n_pt_in_vol_calls;
923
924 // early fail for piv - see if point inside the root level obb
925 // if its not even in the box dont bother doing anything else
926 ErrorCode rval = point_in_box( volume, xyz, result );
927 if( !result )
928 {
929 result = 0;
930 return MB_SUCCESS;
931 }
932
933 // get OBB Tree for volume
934 EntityHandle root;
935 rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to find the volume's obb tree root" );
936
937 // Don't recreate these every call. These cannot be the same as the ray_fire
938 // vectors because both are used simultaneously.
939 std::vector< double > dists;
940 std::vector< EntityHandle > surfs;
941 std::vector< EntityHandle > facets;
942 std::vector< int > dirs;
943
944 // if uvw is not given or is full of zeros, use a random direction
945 double u = 0, v = 0, w = 0;
946
947 if( uvw )
948 {
949 u = uvw[0];
950 v = uvw[1], w = uvw[2];
951 }
952
953 if( u == 0 && v == 0 && w == 0 )
954 {
955 u = rand();
956 v = rand();
957 w = rand();
958 const double magnitude = sqrt( u * u + v * v + w * w );
959 u /= magnitude;
960 v /= magnitude;
961 w /= magnitude;
962 }
963
964 const double ray_direction[] = { u, v, w };
965
966 // if overlaps, ray must be cast to infinity and all RTIs must be returned
967 const double large = 1e15;
968 double ray_length = large;
969
970 // If overlaps occur, the pt is inside if traveling along the ray from the
971 // origin, there are ever more exits than entrances. In lieu of implementing
972 // that, all intersections to infinity are required if overlaps occur (expensive)
973 int min_tolerance_intersections;
974 if( 0 != overlapThickness )
975 {
976 min_tolerance_intersections = -1;
977 // only the first intersection is needed if overlaps do not occur (cheap)
978 }
979 else
980 {
981 min_tolerance_intersections = 1;
982 }
983
984 // Get intersection(s) of forward and reverse orientation. Do not return
985 // glancing intersections or previous facets.
986 GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), xyz, ray_direction, numericalPrecision,
987 min_tolerance_intersections, &root, &volume, &senseTag, NULL,
988 history ? &( history->prev_facets ) : NULL );
989
990 OrientedBoxTreeTool::IntersectSearchWindow search_win( &ray_length, (double*)NULL );
991 rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, xyz,
992 ray_direction, search_win, int_reg_ctxt );MB_CHK_SET_ERR( rval, "Ray fire query failed" );
993
994 // determine orientation of all intersections
995 // 1 for entering, 0 for leaving, -1 for tangent
996 // Tangent intersections are not returned from ray_tri_intersect.
997 dirs.resize( dists.size() );
998 for( unsigned i = 0; i < dists.size(); ++i )
999 {
1000 rval = boundary_case( volume, dirs[i], u, v, w, facets[i], surfs[i] );MB_CHK_SET_ERR( rval, "Failed to resolve boundary case" );
1001 }
1002
1003 // count all crossings
1004 if( 0 != overlapThickness )
1005 {
1006 int sum = 0;
1007 for( unsigned i = 0; i < dirs.size(); ++i )
1008 {
1009 if( 1 == dirs[i] )
1010 sum += 1; // +1 for entering
1011 else if( 0 == dirs[i] )
1012 sum -= 1; // -1 for leaving
1013 else if( -1 == dirs[i] )
1014 { // 0 for tangent
1015 std::cout << "direction==tangent" << std::endl;
1016 sum += 0;
1017 }
1018 else
1019 {
1020 MB_SET_ERR( MB_FAILURE, "Error: unknown direction" );
1021 }
1022 }
1023
1024 // inside/outside depends on the sum
1025 if( 0 < sum )
1026 result = 0; // pt is outside (for all vols)
1027 else if( 0 > sum )
1028 result = 1; // pt is inside (for all vols)
1029 else if( geomTopoTool->is_implicit_complement( volume ) )
1030 result = 1; // pt is inside (for impl_compl_vol)
1031 else
1032 result = 0; // pt is outside (for all other vols)
1033
1034 // Only use the first crossing
1035 }
1036 else
1037 {
1038 if( dirs.empty() )
1039 {
1040 result = 0; // pt is outside
1041 }
1042 else
1043 {
1044 int smallest = std::min_element( dists.begin(), dists.end() ) - dists.begin();
1045 if( 1 == dirs[smallest] )
1046 result = 0; // pt is outside
1047 else if( 0 == dirs[smallest] )
1048 result = 1; // pt is inside
1049 else if( -1 == dirs[smallest] )
1050 {
1051 // Should not be here because Plucker ray-triangle test does not
1052 // return coplanar rays as intersections.
1053 std::cerr << "direction==tangent" << std::endl;
1054 result = -1;
1055 }
1056 else
1057 {
1058 MB_SET_ERR( MB_FAILURE, "Error: unknown direction" );
1059 }
1060 }
1061 }
1062
1063 #ifdef MB_GQT_DEBUG
1064 std::cout << "pt_in_vol: result=" << result << " xyz=" << xyz[0] << " " << xyz[1] << " " << xyz[2] << " uvw=" << u
1065 << " " << v << " " << w << " vol_id=" << volume
1066 << std::endl; // todo: use geomtopotool to get id by entity handle
1067 #endif
1068
1069 return MB_SUCCESS;
1070 }
References boundary_case(), counting, ErrorCode, geomTopoTool, moab::GeomTopoTool::get_root(), moab::GeomTopoTool::is_implicit_complement(), MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, n_pt_in_vol_calls, numericalPrecision, moab::GeomTopoTool::obb_tree(), overlapThickness, point_in_box(), moab::GeomQueryTool::RayHistory::prev_facets, moab::OrientedBoxTreeTool::ray_intersect_sets(), senseTag, and moab::sum().
Referenced by moab::GeomTopoTool::A_is_in_B(), find_volume_slow(), and ray_fire().
| ErrorCode moab::GeomQueryTool::point_in_volume_slow | ( | const EntityHandle | volume, |
| const double | xyz[3], | ||
| int & | result | ||
| ) |
Robust test if a point is inside or outside a volume using unit sphere area method.
This test may be more robust that the standard point_in_volume, but is much slower. It does not detect 'on boundary' situations as point_in_volume does.
| volume | The volume to test |
| xyz | The location to test for volume containment |
| result | Set to 0 if xyz it outside volume, 1 if inside. |
Definition at line 1141 of file GeomQueryTool.cpp.
1142 {
1143 ErrorCode rval;
1144 Range faces;
1145 std::vector< EntityHandle > surfs;
1146 std::vector< int > senses;
1147 double sum = 0.0;
1148 const CartVect point( xyz );
1149
1150 rval = MBI->get_child_meshsets( volume, surfs );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" );
1151
1152 senses.resize( surfs.size() );
1153 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" );
1154
1155 for( unsigned i = 0; i < surfs.size(); ++i )
1156 {
1157 if( !senses[i] ) // skip non-manifold surfaces
1158 continue;
1159
1160 double surf_area = 0.0, face_area;
1161 faces.clear();
1162 rval = MBI->get_entities_by_dimension( surfs[i], 2, faces );MB_CHK_SET_ERR( rval, "Failed to get the surface entities by dimension" );
1163
1164 for( Range::iterator j = faces.begin(); j != faces.end(); ++j )
1165 {
1166 rval = poly_solid_angle( *j, point, face_area );MB_CHK_SET_ERR( rval, "Failed to determin the polygon's solid angle" );
1167
1168 surf_area += face_area;
1169 }
1170
1171 sum += senses[i] * surf_area;
1172 }
1173
1174 result = fabs( sum ) > 2.0 * M_PI;
1175 return MB_SUCCESS;
1176 }
References moab::Range::begin(), moab::Range::clear(), moab::Range::end(), ErrorCode, geomTopoTool, moab::Interface::get_child_meshsets(), moab::Interface::get_entities_by_dimension(), moab::GeomTopoTool::get_surface_senses(), MB_CHK_SET_ERR, MB_SUCCESS, MBI, poly_solid_angle(), and moab::sum().
|
private |
get the solid angle projected by a facet on a unit sphere around a point
Definition at line 1560 of file GeomQueryTool.cpp.
1561 {
1562 ErrorCode rval;
1563
1564 // Get connectivity
1565 const EntityHandle* conn;
1566 int len;
1567 rval = MBI->get_connectivity( face, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the polygon" );
1568
1569 // Allocate space to store vertices
1570 CartVect coords_static[4];
1571 std::vector< CartVect > coords_dynamic;
1572 CartVect* coords = coords_static;
1573 if( (unsigned)len > ( sizeof( coords_static ) / sizeof( coords_static[0] ) ) )
1574 {
1575 coords_dynamic.resize( len );
1576 coords = &coords_dynamic[0];
1577 }
1578
1579 // get coordinates
1580 rval = MBI->get_coords( conn, len, coords->array() );MB_CHK_SET_ERR( rval, "Failed to get the coordinates of the polygon vertices" );
1581
1582 // calculate normal
1583 CartVect norm( 0.0 ), v1, v0 = coords[1] - coords[0];
1584 for( int i = 2; i < len; ++i )
1585 {
1586 v1 = coords[i] - coords[0];
1587 norm += v0 * v1;
1588 v0 = v1;
1589 }
1590
1591 // calculate area
1592 double s, ang;
1593 area = 0.0;
1594 CartVect r, n1, n2, b, a = coords[len - 1] - coords[0];
1595 for( int i = 0; i < len; ++i )
1596 {
1597 r = coords[i] - point;
1598 b = coords[( i + 1 ) % len] - coords[i];
1599 n1 = a * r; // = norm1 (magnitude is important)
1600 n2 = r * b; // = norm2 (magnitude is important)
1601 s = ( n1 % n2 ) / ( n1.length() * n2.length() ); // = cos(angle between norm1,norm2)
1602 ang = s <= -1.0 ? M_PI : s >= 1.0 ? 0.0 : acos( s ); // = acos(s)
1603 s = ( b * a ) % norm; // =orientation of triangle wrt point
1604 area += s > 0.0 ? M_PI - ang : M_PI + ang;
1605 a = -b;
1606 }
1607
1608 area -= M_PI * ( len - 2 );
1609 if( ( norm % r ) > 0 ) area = -area;
1610 return MB_SUCCESS;
1611 }
References moab::CartVect::array(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::CartVect::length(), MB_CHK_SET_ERR, MB_SUCCESS, and MBI.
Referenced by point_in_volume_slow().
| ErrorCode moab::GeomQueryTool::ray_fire | ( | const EntityHandle | volume, |
| const double | ray_start[3], | ||
| const double | ray_dir[3], | ||
| EntityHandle & | next_surf, | ||
| double & | next_surf_dist, | ||
| RayHistory * | history = NULL, |
||
| double | dist_limit = 0, |
||
| int | ray_orientation = 1, |
||
| OrientedBoxTreeTool::TrvStats * | stats = NULL |
||
| ) |
find the next surface crossing from a given point in a given direction
This is the primary method to enable ray tracing through a geometry. Given a volume and a ray, it determines the distance to the nearest intersection with a bounding surface of that volume and returns that distance and the EntityHandle indicating on which surface that intersection occurs. The caller can compute the location of the intersection by adding the distance to the ray.
When a series of calls to this function are made along the same ray (e.g. for the purpose of tracking a ray through several volumes), the optional history argument should be given. The history prevents previously intersected facets from being intersected again. A single history should be used as long as a ray is proceeding forward without changing direction. This situation is sometimes referred to as "streaming."
If a ray changes direction at an intersection site, the caller should call reset_to_last_intersection() on the history object before the next ray fire.
| volume | The volume at which to fire the ray |
| ray_start | An array of x,y,z coordinates from which to start the ray. |
| ray_dir | An array of x,y,z coordinates indicating the direction of the ray. Must be of unit length. |
| next_surf | Output parameter indicating the next surface intersected by the ray. If no intersection is found, will be set to 0. |
| next_surf_dist | Output parameter indicating distance to next_surf. If next_surf is 0, this value is undefined and should not be used. |
| history | Optional RayHistory object. If provided, the facets in the history are assumed to not intersect with the given ray. The facet intersected by this query will also be added to the history. |
| dist_limit | Optional distance limit. If provided and > 0, no intersections at a distance further than this value will be returned. |
| ray_orientation | Optional ray orientation. If provided determines intersections along the normal provided, e.g. if -1 allows intersections back along the the ray direction, Default is 1, i.e. exit intersections |
| stats | Optional TrvStats object used to measure performance of underlying OBB ray-firing query. See OrientedBoxTreeTool.hpp for details. |
Definition at line 738 of file GeomQueryTool.cpp.
747 {
748
749 // take some stats that are independent of nps
750 if( counting )
751 {
752 ++n_ray_fire_calls;
753 if( 0 == n_ray_fire_calls % 10000000 )
754 {
755 std::cout << "n_ray_fires=" << n_ray_fire_calls << " n_pt_in_vols=" << n_pt_in_vol_calls << std::endl;
756 }
757 }
758
759 #ifdef MB_GQT_DEBUG
760 std::cout << "ray_fire:"
761 << " xyz=" << point[0] << " " << point[1] << " " << point[2] << " uvw=" << dir[0] << " " << dir[1] << " "
762 << dir[2] << " entity_handle=" << volume << std::endl;
763 #endif
764
765 const double huge_val = std::numeric_limits< double >::max();
766 double dist_limit = huge_val;
767 if( user_dist_limit > 0 ) dist_limit = user_dist_limit;
768
769 // don't recreate these every call
770 std::vector< double > dists;
771 std::vector< EntityHandle > surfs;
772 std::vector< EntityHandle > facets;
773
774 EntityHandle root;
775 ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the obb tree root of the volume" );
776
777 // check behind the ray origin for intersections
778 double neg_ray_len;
779 if( 0 == overlapThickness )
780 {
781 neg_ray_len = -numericalPrecision;
782 }
783 else
784 {
785 neg_ray_len = -overlapThickness;
786 }
787
788 // optionally, limit the nonneg_ray_len with the distance to next collision.
789 double nonneg_ray_len = dist_limit;
790
791 // the nonneg_ray_len should not be less than -neg_ray_len, or an overlap
792 // may be missed due to optimization within ray_intersect_sets
793 if( nonneg_ray_len < -neg_ray_len ) nonneg_ray_len = -neg_ray_len;
794 if( 0 > nonneg_ray_len || 0 <= neg_ray_len )
795 {
796 MB_SET_ERR( MB_FAILURE, "Incorrect ray length provided" );
797 }
798
799 // min_tolerance_intersections is passed but not used in this call
800 const int min_tolerance_intersections = 0;
801
802 // numericalPrecision is used for box.intersect_ray and find triangles in the
803 // neighborhood of edge/node intersections.
804 GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), point, dir, numericalPrecision, min_tolerance_intersections,
805 &root, &volume, &senseTag, &ray_orientation,
806 history ? &( history->prev_facets ) : NULL );
807
808 OrientedBoxTreeTool::IntersectSearchWindow search_win( &nonneg_ray_len, &neg_ray_len );
809 rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, point, dir,
810 search_win, int_reg_ctxt, stats );
811
812 MB_CHK_SET_ERR( rval, "Ray query failed" );
813
814 // If no distances are returned, the particle is lost unless the physics limit
815 // is being used. If the physics limit is being used, there is no way to tell
816 // if the particle is lost. To avoid ambiguity, DO NOT use the distance limit
817 // unless you know lost particles do not occur.
818 if( dists.empty() )
819 {
820 next_surf = 0;
821 #ifdef MB_GQT_DEBUG
822 std::cout << " next_surf=0 dist=(undef)" << std::endl;
823 #endif
824 return MB_SUCCESS;
825 }
826
827 // Assume that a (neg, nonneg) pair of RTIs could be returned,
828 // however, only one or the other may exist. dists[] may be populated, but
829 // intersections are ONLY indicated by nonzero surfs[] and facets[].
830 if( 2 != dists.size() || 2 != facets.size() )
831 {
832 MB_SET_ERR( MB_FAILURE, "Incorrect number of facets/distances" );
833 }
834 if( 0.0 < dists[0] || 0.0 > dists[1] )
835 {
836 MB_SET_ERR( MB_FAILURE, "Invalid intersection distance signs" );
837 }
838
839 // If both negative and nonnegative RTIs are returned, the negative RTI must
840 // closer to the origin.
841 if( ( 0 != facets[0] && 0 != facets[1] ) && ( -dists[0] > dists[1] ) )
842 {
843 MB_SET_ERR( MB_FAILURE, "Invalid intersection distance values" );
844 }
845
846 // If an RTI is found at negative distance, perform a PMT to see if the
847 // particle is inside an overlap.
848 int exit_idx = -1;
849 if( 0 != facets[0] )
850 {
851 // get the next volume
852 std::vector< EntityHandle > vols;
853 EntityHandle nx_vol;
854 rval = MBI->get_parent_meshsets( surfs[0], vols );MB_CHK_SET_ERR( rval, "Failed to get the parent meshsets" );
855 if( 2 != vols.size() )
856 {
857 MB_SET_ERR( MB_FAILURE, "Invaid number of parent volumes found" );
858 }
859 if( vols.front() == volume )
860 {
861 nx_vol = vols.back();
862 }
863 else
864 {
865 nx_vol = vols.front();
866 }
867 // Check to see if the point is actually in the next volume.
868 // The list of previous facets is used to topologically identify the
869 // "on_boundary" result of the PMT. This avoids a test that uses proximity
870 // (a tolerance).
871 int result;
872 rval = point_in_volume( nx_vol, point, result, dir, history );MB_CHK_SET_ERR( rval, "Point in volume query failed" );
873 if( 1 == result ) exit_idx = 0;
874 }
875
876 // if the negative distance is not the exit, try the nonnegative distance
877 if( -1 == exit_idx && 0 != facets[1] ) exit_idx = 1;
878
879 // if the exit index is still unknown, the particle is lost
880 if( -1 == exit_idx )
881 {
882 next_surf = 0;
883 #ifdef MB_GQT_DEBUG
884 std::cout << "next surf hit = 0, dist = (undef)" << std::endl;
885 #endif
886 return MB_SUCCESS;
887 }
888
889 // return the intersection
890 next_surf = surfs[exit_idx];
891 next_surf_dist = ( 0 > dists[exit_idx] ? 0 : dists[exit_idx] );
892
893 if( history )
894 {
895 history->prev_facets.push_back( facets[exit_idx] );
896 }
897
898 #ifdef MB_GQT_DEBUG
899 if( 0 > dists[exit_idx] )
900 {
901 std::cout << " OVERLAP track length=" << dists[exit_idx] << std::endl;
902 }
903 std::cout << " next_surf = " << next_surf // todo: use geomtopotool to get id by entity handle
904 << ", dist = " << next_surf_dist << " new_pt=";
905 for( int i = 0; i < 3; ++i )
906 {
907 std::cout << point[i] + dir[i] * next_surf_dist << " ";
908 }
909 std::cout << std::endl;
910 #endif
911
912 return MB_SUCCESS;
913 }
References counting, ErrorCode, geomTopoTool, moab::Interface::get_parent_meshsets(), moab::GeomTopoTool::get_root(), MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, MBI, n_pt_in_vol_calls, n_ray_fire_calls, numericalPrecision, moab::GeomTopoTool::obb_tree(), overlapThickness, point_in_volume(), moab::GeomQueryTool::RayHistory::prev_facets, moab::OrientedBoxTreeTool::ray_intersect_sets(), and senseTag.
| void moab::GeomQueryTool::set_numerical_precision | ( | double | new_precision | ) |
Attempt to set a new numerical precision , first checking for sanity Use of this function is discouraged.
Definition at line 1627 of file GeomQueryTool.cpp.
1628 {
1629
1630 if( new_precision <= 0 || new_precision > 1 )
1631 {
1632 std::cerr << "Invalid numerical_precision = " << numericalPrecision << std::endl;
1633 }
1634 else
1635 {
1636 numericalPrecision = new_precision;
1637 }
1638 if( verbose ) std::cout << "Set numerical precision = " << numericalPrecision << std::endl;
1639 }
References numericalPrecision, and verbose.
| void moab::GeomQueryTool::set_overlap_thickness | ( | double | new_overlap_thickness | ) |
State object used in calls to ray_fire()
Storage for the "history" of a ray. This represents the surface facets that the ray is known to have crossed, which cannot be crossed again as long as the ray does not change direction. It is intended to be used with a series of consecutive calls to ray_fire(), in which a ray passes over potentially many surfaces. Attempt to set a new overlap thickness tolerance, first checking for sanity
Definition at line 1613 of file GeomQueryTool.cpp.
1614 {
1615
1616 if( new_thickness < 0 || new_thickness > 100 )
1617 {
1618 std::cerr << "Invalid overlap_thickness = " << new_thickness << std::endl;
1619 }
1620 else
1621 {
1622 overlapThickness = new_thickness;
1623 }
1624 if( verbose ) std::cout << "Set overlap thickness = " << overlapThickness << std::endl;
1625 }
References overlapThickness, and verbose.
|
inline |
Definition at line 366 of file GeomQueryTool.hpp.
367 { 368 verbose = value; 369 }
References verbose.
| ErrorCode moab::GeomQueryTool::test_volume_boundary | ( | const EntityHandle | volume, |
| const EntityHandle | surface, | ||
| const double | xyz[3], | ||
| const double | uvw[3], | ||
| int & | result, | ||
| const RayHistory * | history = NULL |
||
| ) |
Given a ray starting at a surface of a volume, check whether the ray enters or exits the volume.
This function is most useful for rays that change directions at a surface crossing. It can be used to check whether a direction change redirects the ray back into the originating volume.
| volume | The volume to test |
| surface | A surface on volume |
| xyz | A point location on surface |
| uvw | A (unit) direction vector |
| result | Set to 1 if ray is entering volume, or 0 if it is leaving |
| history | Optional ray history object from a previous call to ray_fire. If present and non-empty, the history is used to look up the surface facet at which the ray begins. Absent a history, the facet nearest to xyz will be looked up. The history should always be provided if available, as it avoids the computational expense of a nearest-facet query. |
Definition at line 1103 of file GeomQueryTool.cpp.
1109 {
1110 ErrorCode rval;
1111 int dir;
1112
1113 if( history && history->prev_facets.size() )
1114 {
1115 // the current facet is already available
1116 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" );
1117 }
1118 else
1119 {
1120 // look up nearest facet
1121
1122 // Get OBB Tree for surface
1123 EntityHandle root;
1124 rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's OBB tree root" );
1125
1126 // Get closest triangle on surface
1127 const CartVect point( xyz );
1128 CartVect nearest;
1129 EntityHandle facet_out;
1130 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" );
1131
1132 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" );
1133 }
1134
1135 result = dir;
1136
1137 return MB_SUCCESS;
1138 }
References moab::CartVect::array(), boundary_case(), moab::OrientedBoxTreeTool::closest_to_location(), ErrorCode, geomTopoTool, moab::GeomTopoTool::get_root(), MB_CHK_SET_ERR, MB_SUCCESS, moab::GeomTopoTool::obb_tree(), and moab::GeomQueryTool::RayHistory::prev_facets.
|
private |
Definition at line 402 of file GeomQueryTool.hpp.
Referenced by GeomQueryTool(), point_in_volume(), and ray_fire().
|
private |
Definition at line 397 of file GeomQueryTool.hpp.
Referenced by boundary_case(), closest_to_location(), find_volume(), find_volume_slow(), GeomQueryTool(), get_normal(), gttool(), initialize(), measure_volume(), point_in_box(), point_in_volume(), point_in_volume_slow(), ray_fire(), test_volume_boundary(), and ~GeomQueryTool().
|
private |
Definition at line 400 of file GeomQueryTool.hpp.
Referenced by boundary_case(), find_volume(), GeomQueryTool(), get_normal(), measure_area(), measure_volume(), moab_instance(), point_in_volume_slow(), poly_solid_angle(), and ray_fire().
|
private |
Definition at line 403 of file GeomQueryTool.hpp.
Referenced by GeomQueryTool(), point_in_volume(), and ray_fire().
|
private |
Definition at line 404 of file GeomQueryTool.hpp.
Referenced by GeomQueryTool(), and ray_fire().
|
private |
Definition at line 405 of file GeomQueryTool.hpp.
Referenced by find_volume(), GeomQueryTool(), get_normal(), get_numerical_precision(), point_in_volume(), ray_fire(), and set_numerical_precision().
|
private |
Definition at line 401 of file GeomQueryTool.hpp.
Referenced by GeomQueryTool().
|
private |
Definition at line 405 of file GeomQueryTool.hpp.
Referenced by GeomQueryTool(), get_overlap_thickness(), point_in_volume(), ray_fire(), and set_overlap_thickness().
|
private |
Definition at line 399 of file GeomQueryTool.hpp.
Referenced by ~GeomQueryTool().
|
private |
Definition at line 406 of file GeomQueryTool.hpp.
Referenced by GeomQueryTool(), point_in_volume(), and ray_fire().
|
private |
Definition at line 398 of file GeomQueryTool.hpp.
Referenced by get_verbosity(), set_numerical_precision(), set_overlap_thickness(), and set_verbosity().
1.9.1.