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

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.

 : verbose( false ), owns_gtt( true )
{
 geomTopoTool = new GeomTopoTool( impl, find_geomsets, modelRootSet, p_rootSets_vector, restore_rootSets );
 
 senseTag = geomTopoTool->get_sense_tag();
 
 obbTreeTool = geomTopoTool->obb_tree();
 MBI = geomTopoTool->get_moab_instance();
 
 counting = trace_counting;
 overlapThickness = overlap_thickness;
 numericalPrecision = numerical_precision;
 
 // reset query counters
 n_pt_in_vol_calls = 0;
 n_ray_fire_calls = 0;
}

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.

 : verbose( false ), owns_gtt( false )
{
 
 geomTopoTool = geomtopotool;
 
 senseTag = geomTopoTool->get_sense_tag();
 
 obbTreeTool = geomTopoTool->obb_tree();
 MBI = geomTopoTool->get_moab_instance();
 
 counting = trace_counting;
 overlapThickness = overlap_thickness;
 numericalPrecision = numerical_precision;
 
 // reset query counters
 n_pt_in_vol_calls = 0;
 n_ray_fire_calls = 0;
}

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.

{
 if( owns_gtt )
 {
 delete geomTopoTool;
 }
}

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.

{
 ErrorCode rval;
 
 // test to see if uvw is provided
 if( u <= 1.0 && v <= 1.0 && w <= 1.0 )
 {
 
 const CartVect ray_vector( u, v, w );
 CartVect coords[3], normal( 0.0 );
 const EntityHandle* conn;
 int len, sense_out;
 
 rval = MBI->get_connectivity( facet, conn, len );MB_CHK_SET_ERR( rval, "Failed to get the triangle's connectivity" );
 if( 3 != len )
 {
 MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
 }
 
 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" );
 
 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" );
 
 coords[1] -= coords[0];
 coords[2] -= coords[0];
 normal = sense_out * ( coords[1] * coords[2] );
 
 double sense = ray_vector % normal;
 
 if( sense < 0.0 )
 {
 result = 1; // inside or entering
 }
 else if( sense > 0.0 )
 {
 result = 0; // outside or leaving
 }
 else if( sense == 0.0 )
 {
 result = -1; // tangent, therefore on boundary
 }
 else
 {
 result = -1; // failure
 MB_SET_ERR( MB_FAILURE, "Failed to resolve boundary case" );
 }
 
 // if uvw not provided, return on_boundary.
 }
 else
 {
 result = -1; // on boundary
 return MB_SUCCESS;
 }
 
 return MB_SUCCESS;
}

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

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.

{
 // Get OBB Tree for volume
 EntityHandle root;
 ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's obb tree root" );
 
 // Get closest triangles in volume
 const CartVect point( coords );
 CartVect nearest;
 EntityHandle facet_out;
 
 rval = geomTopoTool->obb_tree()->closest_to_location( point.array(), root, nearest.array(), facet_out,
 closest_surface );MB_CHK_SET_ERR( rval, "Failed to get the closest intersection to location" );
 // calculate distance between point and nearest facet
 result = ( point - nearest ).length();
 
 return MB_SUCCESS;
}

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

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.

{
 ErrorCode rval;
 volume = 0;
 
 EntityHandle global_surf_tree_root = geomTopoTool->get_one_vol_root();
 
 // fast check - make sure point is in the implicit complement bounding box
 int ic_result;
 EntityHandle ic_handle;
 rval = geomTopoTool->get_implicit_complement( ic_handle );MB_CHK_SET_ERR( rval, "Failed to get the implicit complement handle" );
 
 rval = point_in_box( ic_handle, xyz, ic_result );MB_CHK_SET_ERR( rval, "Failed to check implicit complement for containment" );
 if( ic_result == 0 )
 {
 volume = 0;
 return MB_ENTITY_NOT_FOUND;
 }
 
 // if geomTopoTool doesn't have a global tree, use a loop over vols (slow)
 if( !global_surf_tree_root )
 {
 rval = find_volume_slow( xyz, volume, dir );
 return rval;
 }
 
 moab::CartVect uvw( 0.0 );
 
 if( dir )
 {
 uvw[0] = dir[0];
 uvw[1] = dir[1];
 uvw[2] = dir[2];
 }
 
 if( uvw == 0.0 )
 {
 uvw[0] = rand();
 uvw[1] = rand();
 uvw[2] = rand();
 }
 
 // always normalize direction
 uvw.normalize();
 
 // fire a ray along dir and get surface
 const double huge_val = std::numeric_limits< double >::max();
 double pos_ray_len = huge_val;
 double neg_ray_len = -huge_val;
 
 // RIS output data
 std::vector< double > dists;
 std::vector< EntityHandle > surfs;
 std::vector< EntityHandle > facets;
 
 FindVolumeIntRegCtxt find_vol_reg_ctxt;
 OrientedBoxTreeTool::IntersectSearchWindow search_win( &pos_ray_len, &neg_ray_len );
 rval =
 geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, global_surf_tree_root, numericalPrecision,
 xyz, uvw.array(), search_win, find_vol_reg_ctxt );MB_CHK_SET_ERR( rval, "Failed in global tree ray fire" );
 
 // if there was no intersection, no volume is found
 if( surfs.size() == 0 || surfs[0] == 0 )
 {
 volume = 0;
 return MB_ENTITY_NOT_FOUND;
 }
 
 // get the positive distance facet, surface hit
 EntityHandle facet = facets[0];
 EntityHandle surf = surfs[0];
 
 // get these now, we're going to use them no matter what
 EntityHandle fwd_vol, bwd_vol;
 rval = geomTopoTool->get_surface_senses( surf, fwd_vol, bwd_vol );MB_CHK_SET_ERR( rval, "Failed to get sense data" );
 EntityHandle parent_vols[2];
 parent_vols[0] = fwd_vol;
 parent_vols[1] = bwd_vol;
 
 // get triangle normal
 std::vector< EntityHandle > conn;
 CartVect coords[3];
 rval = MBI->get_connectivity( &facet, 1, conn );MB_CHK_SET_ERR( rval, "Failed to get triangle connectivity" );
 
 rval = MBI->get_coords( &conn[0], 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get triangle coordinates" );
 
 CartVect normal = ( coords[1] - coords[0] ) * ( coords[2] - coords[0] );
 normal.normalize();
 
 // reverse direction if a hit in the negative direction is found
 if( dists[0] < 0 )
 {
 uvw *= -1;
 }
 
 // if this is a "forward" intersection return the first sense entity
 // otherwise return the second, "reverse" sense entity
 double dot_prod = uvw % normal;
 int idx = dot_prod > 0.0 ? 0 : 1;
 
 if( dot_prod == 0.0 )
 {
 std::cerr << "Tangent dot product in find_volume. Shouldn't be here." << std::endl;
 volume = 0;
 return MB_FAILURE;
 }
 
 volume = parent_vols[idx];
 
 return MB_SUCCESS;
}

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

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.

{
 ErrorCode rval;
 volume = 0;
 // get all volumes
 Range all_vols;
 rval = geomTopoTool->get_gsets_by_dimension( 3, all_vols );MB_CHK_SET_ERR( rval, "Failed to get all volumes in the model" );
 
 Range::iterator it;
 int result = 0;
 for( it = all_vols.begin(); it != all_vols.end(); it++ )
 {
 rval = point_in_volume( *it, xyz, result, dir );MB_CHK_SET_ERR( rval, "Failed in point in volume loop" );
 if( result )
 {
 volume = *it;
 break;
 }
 }
 return volume ? MB_SUCCESS : MB_ENTITY_NOT_FOUND;
}

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

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.

{
 EntityHandle root;
 ErrorCode rval = geomTopoTool->get_root( surf, root );MB_CHK_SET_ERR( rval, "Failed to get the surface's obb tree root" );
 
 std::vector< EntityHandle > facets;
 
 // if no history or history empty, use nearby facets
 if( !history || ( history->prev_facets.size() == 0 ) )
 {
 rval = geomTopoTool->obb_tree()->closest_to_location( in_pt, root, numericalPrecision, facets );MB_CHK_SET_ERR( rval, "Failed to get closest intersection to location" );
 }
 // otherwise use most recent facet in history
 else
 {
 facets.push_back( history->prev_facets.back() );
 }
 
 CartVect coords[3], normal( 0.0 );
 const EntityHandle* conn;
 int len;
 for( unsigned i = 0; i < facets.size(); ++i )
 {
 rval = MBI->get_connectivity( facets[i], conn, len );MB_CHK_SET_ERR( rval, "Failed to get facet connectivity" );
 if( 3 != len )
 {
 MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
 }
 
 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" );
 
 coords[1] -= coords[0];
 coords[2] -= coords[0];
 normal += coords[1] * coords[2];
 }
 
 normal.normalize();
 normal.get( angle );
 
 return MB_SUCCESS;
}

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.

 {
 return numericalPrecision;
 }

References numericalPrecision.

get_overlap_thickness()

double moab::GeomQueryTool::get_overlap_thickness ( )

inline

Definition at line 381 of file GeomQueryTool.hpp.

 {
 return overlapThickness;
 }

References overlapThickness.

get_verbosity()

bool moab::GeomQueryTool::get_verbosity ( ) const

inline

Definition at line 371 of file GeomQueryTool.hpp.

 {
 return verbose;
 }

References verbose.

gttool()

GeomTopoTool* moab::GeomQueryTool::gttool ( )

inline

Definition at line 386 of file GeomQueryTool.hpp.

 {
 return geomTopoTool;
 }

References geomTopoTool.

initialize()

ErrorCode moab::GeomQueryTool::initialize

(

)

Definition at line 681 of file GeomQueryTool.cpp.

{
 
 ErrorCode rval;
 
 rval = geomTopoTool->find_geomsets();MB_CHK_SET_ERR( rval, "Failed to find geometry sets" );
 
 rval = geomTopoTool->setup_implicit_complement();MB_CHK_SET_ERR( rval, "Couldn't setup the implicit complement" );
 
 rval = geomTopoTool->construct_obb_trees();MB_CHK_SET_ERR( rval, "Failed to construct OBB trees" );
 
 return MB_SUCCESS;
}

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.

{
 // get triangles in surface
 Range triangles;
 ErrorCode rval = MBI->get_entities_by_dimension( surface, 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface entities" );
 if( !triangles.all_of_type( MBTRI ) )
 {
 std::cout << "WARNING: Surface " << surface // todo: use geomtopotool to get id by entity handle
 << " contains non-triangle elements. Area calculation may be incorrect." << std::endl;
 triangles.clear();
 rval = MBI->get_entities_by_type( surface, MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to the surface's triangle entities" );
 }
 
 // calculate sum of area of triangles
 result = 0.0;
 const EntityHandle* conn;
 int len;
 CartVect coords[3];
 for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j )
 {
 rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's connectivity" );
 if( 3 != len )
 {
 MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
 }
 rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's vertex coordinates" );
 
 // calculated area using cross product of triangle edges
 CartVect v1 = coords[1] - coords[0];
 CartVect v2 = coords[2] - coords[0];
 CartVect xp = v1 * v2;
 result += xp.length();
 }
 result *= 0.5;
 return MB_SUCCESS;
}

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.

{
 ErrorCode rval;
 std::vector< EntityHandle > surfaces;
 result = 0.0;
 
 // don't try to calculate volume of implicit complement
 if( geomTopoTool->is_implicit_complement( volume ) )
 {
 result = 1.0;
 return MB_SUCCESS;
 }
 
 // get surfaces from volume
 rval = MBI->get_child_meshsets( volume, surfaces );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" );
 
 // get surface senses
 std::vector< int > senses( surfaces.size() );
 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" );
 
 for( unsigned i = 0; i < surfaces.size(); ++i )
 {
 // skip non-manifold surfaces
 if( !senses[i] ) continue;
 
 // get triangles in surface
 Range triangles;
 rval = MBI->get_entities_by_dimension( surfaces[i], 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" );
 
 if( !triangles.all_of_type( MBTRI ) )
 {
 std::cout << "WARNING: Surface " << surfaces[i] // todo: use geomtopotool to get id by entity handle
 << " contains non-triangle elements. Volume calculation may be incorrect." << std::endl;
 triangles.clear();
 rval = MBI->get_entities_by_type( surfaces[i], MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" );
 }
 
 // calculate signed volume beneath surface (x 6.0)
 double surf_sum = 0.0;
 const EntityHandle* conn;
 int len;
 CartVect coords[3];
 for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j )
 {
 rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the current triangle" );
 if( 3 != len )
 {
 MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
 }
 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" );
 
 coords[1] -= coords[0];
 coords[2] -= coords[0];
 surf_sum += ( coords[0] % ( coords[1] * coords[2] ) );
 }
 result += senses[i] * surf_sum;
 }
 
 result /= 6.0;
 return MB_SUCCESS;
}

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.

Referenced by main(), and summarize_cell_volume_change().

moab_instance()

Interface* moab::GeomQueryTool::moab_instance ( )

inline

Definition at line 391 of file GeomQueryTool.hpp.

 {
 return MBI;
 }

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

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.

{
 double minpt[3];
 double maxpt[3];
 ErrorCode rval = geomTopoTool->get_bounding_coords( volume, minpt, maxpt );MB_CHK_SET_ERR( rval, "Failed to get the bounding coordinates of the volume" );
 
 // early exits
 if( point[0] > maxpt[0] || point[0] < minpt[0] )
 {
 inside = 0;
 return rval;
 }
 if( point[1] > maxpt[1] || point[1] < minpt[1] )
 {
 inside = 0;
 return rval;
 }
 if( point[2] > maxpt[2] || point[2] < minpt[2] )
 {
 inside = 0;
 return rval;
 }
 inside = 1;
 return rval;
}

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

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.

{
 // take some stats that are independent of nps
 if( counting ) ++n_pt_in_vol_calls;
 
 // early fail for piv - see if point inside the root level obb
 // if its not even in the box dont bother doing anything else
 ErrorCode rval = point_in_box( volume, xyz, result );
 if( !result )
 {
 result = 0;
 return MB_SUCCESS;
 }
 
 // get OBB Tree for volume
 EntityHandle root;
 rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to find the volume's obb tree root" );
 
 // Don't recreate these every call. These cannot be the same as the ray_fire
 // vectors because both are used simultaneously.
 std::vector< double > dists;
 std::vector< EntityHandle > surfs;
 std::vector< EntityHandle > facets;
 std::vector< int > dirs;
 
 // if uvw is not given or is full of zeros, use a random direction
 double u = 0, v = 0, w = 0;
 
 if( uvw )
 {
 u = uvw[0];
 v = uvw[1], w = uvw[2];
 }
 
 if( u == 0 && v == 0 && w == 0 )
 {
 u = rand();
 v = rand();
 w = rand();
 const double magnitude = sqrt( u * u + v * v + w * w );
 u /= magnitude;
 v /= magnitude;
 w /= magnitude;
 }
 
 const double ray_direction[] = { u, v, w };
 
 // if overlaps, ray must be cast to infinity and all RTIs must be returned
 const double large = 1e15;
 double ray_length = large;
 
 // If overlaps occur, the pt is inside if traveling along the ray from the
 // origin, there are ever more exits than entrances. In lieu of implementing
 // that, all intersections to infinity are required if overlaps occur (expensive)
 int min_tolerance_intersections;
 if( 0 != overlapThickness )
 {
 min_tolerance_intersections = -1;
 // only the first intersection is needed if overlaps do not occur (cheap)
 }
 else
 {
 min_tolerance_intersections = 1;
 }
 
 // Get intersection(s) of forward and reverse orientation. Do not return
 // glancing intersections or previous facets.
 GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), xyz, ray_direction, numericalPrecision,
 min_tolerance_intersections, &root, &volume, &senseTag, NULL,
 history ? &( history->prev_facets ) : NULL );
 
 OrientedBoxTreeTool::IntersectSearchWindow search_win( &ray_length, (double*)NULL );
 rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, xyz,
 ray_direction, search_win, int_reg_ctxt );MB_CHK_SET_ERR( rval, "Ray fire query failed" );
 
 // determine orientation of all intersections
 // 1 for entering, 0 for leaving, -1 for tangent
 // Tangent intersections are not returned from ray_tri_intersect.
 dirs.resize( dists.size() );
 for( unsigned i = 0; i < dists.size(); ++i )
 {
 rval = boundary_case( volume, dirs[i], u, v, w, facets[i], surfs[i] );MB_CHK_SET_ERR( rval, "Failed to resolve boundary case" );
 }
 
 // count all crossings
 if( 0 != overlapThickness )
 {
 int sum = 0;
 for( unsigned i = 0; i < dirs.size(); ++i )
 {
 if( 1 == dirs[i] )
 sum += 1; // +1 for entering
 else if( 0 == dirs[i] )
 sum -= 1; // -1 for leaving
 else if( -1 == dirs[i] )
 { // 0 for tangent
 std::cout << "direction==tangent" << std::endl;
 sum += 0;
 }
 else
 {
 MB_SET_ERR( MB_FAILURE, "Error: unknown direction" );
 }
 }
 
 // inside/outside depends on the sum
 if( 0 < sum )
 result = 0; // pt is outside (for all vols)
 else if( 0 > sum )
 result = 1; // pt is inside (for all vols)
 else if( geomTopoTool->is_implicit_complement( volume ) )
 result = 1; // pt is inside (for impl_compl_vol)
 else
 result = 0; // pt is outside (for all other vols)
 
 // Only use the first crossing
 }
 else
 {
 if( dirs.empty() )
 {
 result = 0; // pt is outside
 }
 else
 {
 int smallest = std::min_element( dists.begin(), dists.end() ) - dists.begin();
 if( 1 == dirs[smallest] )
 result = 0; // pt is outside
 else if( 0 == dirs[smallest] )
 result = 1; // pt is inside
 else if( -1 == dirs[smallest] )
 {
 // Should not be here because Plucker ray-triangle test does not
 // return coplanar rays as intersections.
 std::cerr << "direction==tangent" << std::endl;
 result = -1;
 }
 else
 {
 MB_SET_ERR( MB_FAILURE, "Error: unknown direction" );
 }
 }
 }
 
#ifdef MB_GQT_DEBUG
 std::cout << "pt_in_vol: result=" << result << " xyz=" << xyz[0] << " " << xyz[1] << " " << xyz[2] << " uvw=" << u
 << " " << v << " " << w << " vol_id=" << volume
 << std::endl; // todo: use geomtopotool to get id by entity handle
#endif
 
 return MB_SUCCESS;
}

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

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.

{
 ErrorCode rval;
 Range faces;
 std::vector< EntityHandle > surfs;
 std::vector< int > senses;
 double sum = 0.0;
 const CartVect point( xyz );
 
 rval = MBI->get_child_meshsets( volume, surfs );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" );
 
 senses.resize( surfs.size() );
 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" );
 
 for( unsigned i = 0; i < surfs.size(); ++i )
 {
 if( !senses[i] ) // skip non-manifold surfaces
 continue;
 
 double surf_area = 0.0, face_area;
 faces.clear();
 rval = MBI->get_entities_by_dimension( surfs[i], 2, faces );MB_CHK_SET_ERR( rval, "Failed to get the surface entities by dimension" );
 
 for( Range::iterator j = faces.begin(); j != faces.end(); ++j )
 {
 rval = poly_solid_angle( *j, point, face_area );MB_CHK_SET_ERR( rval, "Failed to determin the polygon's solid angle" );
 
 surf_area += face_area;
 }
 
 sum += senses[i] * surf_area;
 }
 
 result = fabs( sum ) > 2.0 * M_PI;
 return MB_SUCCESS;
}

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.

{
 ErrorCode rval;
 
 // Get connectivity
 const EntityHandle* conn;
 int len;
 rval = MBI->get_connectivity( face, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the polygon" );
 
 // Allocate space to store vertices
 CartVect coords_static[4];
 std::vector< CartVect > coords_dynamic;
 CartVect* coords = coords_static;
 if( (unsigned)len > ( sizeof( coords_static ) / sizeof( coords_static[0] ) ) )
 {
 coords_dynamic.resize( len );
 coords = &coords_dynamic[0];
 }
 
 // get coordinates
 rval = MBI->get_coords( conn, len, coords->array() );MB_CHK_SET_ERR( rval, "Failed to get the coordinates of the polygon vertices" );
 
 // calculate normal
 CartVect norm( 0.0 ), v1, v0 = coords[1] - coords[0];
 for( int i = 2; i < len; ++i )
 {
 v1 = coords[i] - coords[0];
 norm += v0 * v1;
 v0 = v1;
 }
 
 // calculate area
 double s, ang;
 area = 0.0;
 CartVect r, n1, n2, b, a = coords[len - 1] - coords[0];
 for( int i = 0; i < len; ++i )
 {
 r = coords[i] - point;
 b = coords[( i + 1 ) % len] - coords[i];
 n1 = a * r; // = norm1 (magnitude is important)
 n2 = r * b; // = norm2 (magnitude is important)
 s = ( n1 % n2 ) / ( n1.length() * n2.length() ); // = cos(angle between norm1,norm2)
 ang = s <= -1.0 ? M_PI : s >= 1.0 ? 0.0 : acos( s ); // = acos(s)
 s = ( b * a ) % norm; // =orientation of triangle wrt point
 area += s > 0.0 ? M_PI - ang : M_PI + ang;
 a = -b;
 }
 
 area -= M_PI * ( len - 2 );
 if( ( norm % r ) > 0 ) area = -area;
 return MB_SUCCESS;
}

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

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.

{
 
 // take some stats that are independent of nps
 if( counting )
 {
 ++n_ray_fire_calls;
 if( 0 == n_ray_fire_calls % 10000000 )
 {
 std::cout << "n_ray_fires=" << n_ray_fire_calls << " n_pt_in_vols=" << n_pt_in_vol_calls << std::endl;
 }
 }
 
#ifdef MB_GQT_DEBUG
 std::cout << "ray_fire:"
 << " xyz=" << point[0] << " " << point[1] << " " << point[2] << " uvw=" << dir[0] << " " << dir[1] << " "
 << dir[2] << " entity_handle=" << volume << std::endl;
#endif
 
 const double huge_val = std::numeric_limits< double >::max();
 double dist_limit = huge_val;
 if( user_dist_limit > 0 ) dist_limit = user_dist_limit;
 
 // don't recreate these every call
 std::vector< double > dists;
 std::vector< EntityHandle > surfs;
 std::vector< EntityHandle > facets;
 
 EntityHandle root;
 ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the obb tree root of the volume" );
 
 // check behind the ray origin for intersections
 double neg_ray_len;
 if( 0 == overlapThickness )
 {
 neg_ray_len = -numericalPrecision;
 }
 else
 {
 neg_ray_len = -overlapThickness;
 }
 
 // optionally, limit the nonneg_ray_len with the distance to next collision.
 double nonneg_ray_len = dist_limit;
 
 // the nonneg_ray_len should not be less than -neg_ray_len, or an overlap
 // may be missed due to optimization within ray_intersect_sets
 if( nonneg_ray_len < -neg_ray_len ) nonneg_ray_len = -neg_ray_len;
 if( 0 > nonneg_ray_len || 0 <= neg_ray_len )
 {
 MB_SET_ERR( MB_FAILURE, "Incorrect ray length provided" );
 }
 
 // min_tolerance_intersections is passed but not used in this call
 const int min_tolerance_intersections = 0;
 
 // numericalPrecision is used for box.intersect_ray and find triangles in the
 // neighborhood of edge/node intersections.
 GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), point, dir, numericalPrecision, min_tolerance_intersections,
 &root, &volume, &senseTag, &ray_orientation,
 history ? &( history->prev_facets ) : NULL );
 
 OrientedBoxTreeTool::IntersectSearchWindow search_win( &nonneg_ray_len, &neg_ray_len );
 rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, point, dir,
 search_win, int_reg_ctxt, stats );
 
 MB_CHK_SET_ERR( rval, "Ray query failed" );
 
 // If no distances are returned, the particle is lost unless the physics limit
 // is being used. If the physics limit is being used, there is no way to tell
 // if the particle is lost. To avoid ambiguity, DO NOT use the distance limit
 // unless you know lost particles do not occur.
 if( dists.empty() )
 {
 next_surf = 0;
#ifdef MB_GQT_DEBUG
 std::cout << " next_surf=0 dist=(undef)" << std::endl;
#endif
 return MB_SUCCESS;
 }
 
 // Assume that a (neg, nonneg) pair of RTIs could be returned,
 // however, only one or the other may exist. dists[] may be populated, but
 // intersections are ONLY indicated by nonzero surfs[] and facets[].
 if( 2 != dists.size() || 2 != facets.size() )
 {
 MB_SET_ERR( MB_FAILURE, "Incorrect number of facets/distances" );
 }
 if( 0.0 < dists[0] || 0.0 > dists[1] )
 {
 MB_SET_ERR( MB_FAILURE, "Invalid intersection distance signs" );
 }
 
 // If both negative and nonnegative RTIs are returned, the negative RTI must
 // closer to the origin.
 if( ( 0 != facets[0] && 0 != facets[1] ) && ( -dists[0] > dists[1] ) )
 {
 MB_SET_ERR( MB_FAILURE, "Invalid intersection distance values" );
 }
 
 // If an RTI is found at negative distance, perform a PMT to see if the
 // particle is inside an overlap.
 int exit_idx = -1;
 if( 0 != facets[0] )
 {
 // get the next volume
 std::vector< EntityHandle > vols;
 EntityHandle nx_vol;
 rval = MBI->get_parent_meshsets( surfs[0], vols );MB_CHK_SET_ERR( rval, "Failed to get the parent meshsets" );
 if( 2 != vols.size() )
 {
 MB_SET_ERR( MB_FAILURE, "Invaid number of parent volumes found" );
 }
 if( vols.front() == volume )
 {
 nx_vol = vols.back();
 }
 else
 {
 nx_vol = vols.front();
 }
 // Check to see if the point is actually in the next volume.
 // The list of previous facets is used to topologically identify the
 // "on_boundary" result of the PMT. This avoids a test that uses proximity
 // (a tolerance).
 int result;
 rval = point_in_volume( nx_vol, point, result, dir, history );MB_CHK_SET_ERR( rval, "Point in volume query failed" );
 if( 1 == result ) exit_idx = 0;
 }
 
 // if the negative distance is not the exit, try the nonnegative distance
 if( -1 == exit_idx && 0 != facets[1] ) exit_idx = 1;
 
 // if the exit index is still unknown, the particle is lost
 if( -1 == exit_idx )
 {
 next_surf = 0;
#ifdef MB_GQT_DEBUG
 std::cout << "next surf hit = 0, dist = (undef)" << std::endl;
#endif
 return MB_SUCCESS;
 }
 
 // return the intersection
 next_surf = surfs[exit_idx];
 next_surf_dist = ( 0 > dists[exit_idx] ? 0 : dists[exit_idx] );
 
 if( history )
 {
 history->prev_facets.push_back( facets[exit_idx] );
 }
 
#ifdef MB_GQT_DEBUG
 if( 0 > dists[exit_idx] )
 {
 std::cout << " OVERLAP track length=" << dists[exit_idx] << std::endl;
 }
 std::cout << " next_surf = " << next_surf // todo: use geomtopotool to get id by entity handle
 << ", dist = " << next_surf_dist << " new_pt=";
 for( int i = 0; i < 3; ++i )
 {
 std::cout << point[i] + dir[i] * next_surf_dist << " ";
 }
 std::cout << std::endl;
#endif
 
 return MB_SUCCESS;
}

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.

{
 
 if( new_precision <= 0 || new_precision > 1 )
 {
 std::cerr << "Invalid numerical_precision = " << numericalPrecision << std::endl;
 }
 else
 {
 numericalPrecision = new_precision;
 }
 if( verbose ) std::cout << "Set numerical precision = " << numericalPrecision << std::endl;
}

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.

{
 
 if( new_thickness < 0 || new_thickness > 100 )
 {
 std::cerr << "Invalid overlap_thickness = " << new_thickness << std::endl;
 }
 else
 {
 overlapThickness = new_thickness;
 }
 if( verbose ) std::cout << "Set overlap thickness = " << overlapThickness << std::endl;
}

References overlapThickness, and verbose.

set_verbosity()

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

inline

Definition at line 366 of file GeomQueryTool.hpp.

 {
 verbose = value;
 }

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

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.

{
 ErrorCode rval;
 int dir;
 
 if( history && history->prev_facets.size() )
 {
 // the current facet is already available
 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" );
 }
 else
 {
 // look up nearest facet
 
 // Get OBB Tree for surface
 EntityHandle root;
 rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's OBB tree root" );
 
 // Get closest triangle on surface
 const CartVect point( xyz );
 CartVect nearest;
 EntityHandle facet_out;
 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" );
 
 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" );
 }
 
 result = dir;
 
 return MB_SUCCESS;
}

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

GeomTopoTool* moab::GeomQueryTool::geomTopoTool

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

MBI

Interface* moab::GeomQueryTool::MBI

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

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

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

obbTreeTool

OrientedBoxTreeTool* moab::GeomQueryTool::obbTreeTool

private

Definition at line 401 of file GeomQueryTool.hpp.

Referenced by GeomQueryTool().

overlapThickness

double moab::GeomQueryTool::overlapThickness

private

Definition at line 405 of file GeomQueryTool.hpp.

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

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

Definition at line 398 of file GeomQueryTool.hpp.

Referenced by get_verbosity(), set_numerical_precision(), set_overlap_thickness(), and set_verbosity().

---

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

• GeomQueryTool.hpp
• GeomQueryTool.cpp