Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
moab::GeomQueryTool Class Reference

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 ()
 
GeomTopoToolgttool ()
 
Interfacemoab_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

GeomTopoToolgeomTopoTool
 
bool verbose
 
bool owns_gtt
 
InterfaceMBI
 
OrientedBoxTreeToolobbTreeTool
 
bool counting
 
long long int n_pt_in_vol_calls
 
long long int n_ray_fire_calls
 
double overlapThickness
 
double numericalPrecision
 
Tag senseTag
 

Detailed Description

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:

  • measures of surface area and volume
  • distance to a bounding surface from a point within a volume
  • test the inclusion of a point within a volume
  • find the angle of a surface at a point

A feature of these queries is that there is some support for overlapping geometries.

Definition at line 45 of file GeomQueryTool.hpp.

Constructor & Destructor Documentation

◆ GeomQueryTool() [1/2]

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 
637 
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.

◆ GeomQueryTool() [2/2]

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 
660 
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.

◆ ~GeomQueryTool()

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.

Member Function Documentation

◆ boundary_case()

ErrorCode moab::GeomQueryTool::boundary_case ( EntityHandle  volume,
int &  result,
double  u,
double  v,
double  w,
EntityHandle  facet,
EntityHandle  surface 
)
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().

◆ closest_to_location()

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.

Parameters
volumeVolume to query
pointCoordinates of test point
resultSet 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().

◆ find_volume()

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.

Parameters
xyzThe location to test
volumeSet to volume handle containing the location if found
dirOptional 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().

◆ find_volume_slow()

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

Parameters
xyzThe location to test
volumeSet to volume handle containing the location if found
dirOptional 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().

◆ get_normal()

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.

Parameters
surfSurface on which to get normal
xyzPoint on surf
angleSet to coordinates of surface normal nearest xyz
historyOptional 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.

◆ get_numerical_precision()

double moab::GeomQueryTool::get_numerical_precision ( )
inline

Definition at line 376 of file GeomQueryTool.hpp.

377  {
378  return numericalPrecision;
379  }

References numericalPrecision.

◆ get_overlap_thickness()

double moab::GeomQueryTool::get_overlap_thickness ( )
inline

Definition at line 381 of file GeomQueryTool.hpp.

382  {
383  return overlapThickness;
384  }

References overlapThickness.

◆ get_verbosity()

bool moab::GeomQueryTool::get_verbosity ( ) const
inline

Definition at line 371 of file GeomQueryTool.hpp.

372  {
373  return verbose;
374  }

References verbose.

◆ gttool()

GeomTopoTool* moab::GeomQueryTool::gttool ( )
inline

Definition at line 386 of file GeomQueryTool.hpp.

387  {
388  return geomTopoTool;
389  }

References geomTopoTool.

◆ initialize()

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().

◆ measure_area()

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.

◆ measure_volume()

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.

◆ moab_instance()

Interface* moab::GeomQueryTool::moab_instance ( )
inline

Definition at line 391 of file GeomQueryTool.hpp.

392  {
393  return MBI;
394  }

References MBI.

◆ point_in_box()

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.

Parameters
volumeThe volume to test
pointThe location to test for bounding box containment
insideset 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().

◆ 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).

Parameters
volumeThe volume to test
xyzThe location to test for volume containment
resultSet to 0 if xyz it outside volume, 1 if inside, and -1 if on boundary.
Optionaldirection 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.
historyOptional 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().

◆ point_in_volume_slow()

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.

Parameters
volumeThe volume to test
xyzThe location to test for volume containment
resultSet 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().

◆ poly_solid_angle()

ErrorCode moab::GeomQueryTool::poly_solid_angle ( EntityHandle  face,
const CartVect point,
double &  area 
)
private

get the solid angle projected by a facet on a unit sphere around a point

  • used by point_in_volume_slow

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().

◆ ray_fire()

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.

Parameters
volumeThe volume at which to fire the ray
ray_startAn array of x,y,z coordinates from which to start the ray.
ray_dirAn array of x,y,z coordinates indicating the direction of the ray. Must be of unit length.
next_surfOutput parameter indicating the next surface intersected by the ray. If no intersection is found, will be set to 0.
next_surf_distOutput parameter indicating distance to next_surf. If next_surf is 0, this value is undefined and should not be used.
historyOptional 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_limitOptional distance limit. If provided and > 0, no intersections at a distance further than this value will be returned.
ray_orientationOptional 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
statsOptional 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  {
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.

◆ set_numerical_precision()

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.

◆ set_overlap_thickness()

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.

◆ set_verbosity()

void moab::GeomQueryTool::set_verbosity ( bool  value)
inline

Definition at line 366 of file GeomQueryTool.hpp.

367  {
368  verbose = value;
369  }

References verbose.

◆ test_volume_boundary()

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.

Parameters
volumeThe volume to test
surfaceA surface on volume
xyzA point location on surface
uvwA (unit) direction vector
resultSet to 1 if ray is entering volume, or 0 if it is leaving
historyOptional 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.

Member Data Documentation

◆ counting

bool moab::GeomQueryTool::counting
private

Definition at line 402 of file GeomQueryTool.hpp.

Referenced by GeomQueryTool(), point_in_volume(), and ray_fire().

◆ geomTopoTool

◆ MBI

◆ n_pt_in_vol_calls

long long int moab::GeomQueryTool::n_pt_in_vol_calls
private

Definition at line 403 of file GeomQueryTool.hpp.

Referenced by GeomQueryTool(), point_in_volume(), and ray_fire().

◆ n_ray_fire_calls

long long int moab::GeomQueryTool::n_ray_fire_calls
private

Definition at line 404 of file GeomQueryTool.hpp.

Referenced by GeomQueryTool(), and ray_fire().

◆ numericalPrecision

double moab::GeomQueryTool::numericalPrecision
private

◆ obbTreeTool

OrientedBoxTreeTool* moab::GeomQueryTool::obbTreeTool
private

Definition at line 401 of file GeomQueryTool.hpp.

Referenced by GeomQueryTool().

◆ overlapThickness

double moab::GeomQueryTool::overlapThickness
private

◆ owns_gtt

bool moab::GeomQueryTool::owns_gtt
private

Definition at line 399 of file GeomQueryTool.hpp.

Referenced by ~GeomQueryTool().

◆ senseTag

Tag moab::GeomQueryTool::senseTag
private

Definition at line 406 of file GeomQueryTool.hpp.

Referenced by GeomQueryTool(), point_in_volume(), and ray_fire().

◆ verbose

bool moab::GeomQueryTool::verbose
private

The documentation for this class was generated from the following files: