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
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.
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.
ErrorCode	find_volume (const double xyz[3], EntityHandle &volume, const double *dir=NULL)
	Find volume for a given location.
ErrorCode	find_volume_slow (const double xyz[3], EntityHandle &volume, const double *dir=NULL)
	Find volume for a given location using loop. (slow)
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.
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.
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.
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()
void	set_numerical_precision (double new_precision)
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
ErrorCode	poly_solid_angle (EntityHandle face, const CartVect &point, double &area)
Private Attributes
GeomTopoTool *	geomTopoTool
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

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 625 of file GeomQueryTool.cpp.

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.

    : 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;
}

moab::GeomQueryTool::GeomQueryTool	(	GeomTopoTool *	geomtopotool,
		bool	trace_counting = `false`,
		double	overlap_thickness = `0.`,
		double	numerical_precision = `0.001`
	)

Definition at line 651 of file GeomQueryTool.cpp.

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.

    : 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;
}

moab::GeomQueryTool::~GeomQueryTool ( )

Definition at line 674 of file GeomQueryTool.cpp.

References geomTopoTool, and owns_gtt.

{
    if( owns_gtt )
    {
        delete geomTopoTool;
    }
}

Member Function Documentation

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 1492 of file GeomQueryTool.cpp.

References ErrorCode, geomTopoTool, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::GeomTopoTool::get_sense(), MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, and MBI.

Referenced by point_in_volume(), and test_volume_boundary().

{
    ErrorCode 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;
}

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 1317 of file GeomQueryTool.cpp.

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

{
    // 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;
}

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 1182 of file GeomQueryTool.cpp.

References moab::CartVect::array(), ErrorCode, find_volume_slow(), geomTopoTool, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::GeomTopoTool::get_implicit_complement(), moab::GeomTopoTool::get_one_vol_root(), moab::GeomTopoTool::get_surface_senses(), MB_CHK_SET_ERR, MB_ENTITY_NOT_FOUND, MB_SUCCESS, MBI, moab::CartVect::normalize(), numericalPrecision, moab::GeomTopoTool::obb_tree(), point_in_box(), and moab::OrientedBoxTreeTool::ray_intersect_sets().

{
    ErrorCode 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;
}

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 1294 of file GeomQueryTool.cpp.

References moab::Range::begin(), moab::Range::end(), ErrorCode, geomTopoTool, moab::GeomTopoTool::get_gsets_by_dimension(), MB_CHK_SET_ERR, MB_ENTITY_NOT_FOUND, MB_SUCCESS, and point_in_volume().

Referenced by find_volume().

{
    ErrorCode 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;
}

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 1440 of file GeomQueryTool.cpp.

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

{
    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;
}

double moab::GeomQueryTool::get_numerical_precision ( ) [inline]

Definition at line 366 of file GeomQueryTool.hpp.

References numericalPrecision.

    {
        return numericalPrecision;
    }

double moab::GeomQueryTool::get_overlap_thickness ( ) [inline]

Definition at line 371 of file GeomQueryTool.hpp.

References overlapThickness.

    {
        return overlapThickness;
    }

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

Definition at line 376 of file GeomQueryTool.hpp.

References geomTopoTool.

    {
        return geomTopoTool;
    }

ErrorCode moab::GeomQueryTool::initialize ( )

Definition at line 682 of file GeomQueryTool.cpp.

References moab::GeomTopoTool::construct_obb_trees(), ErrorCode, moab::GeomTopoTool::find_geomsets(), geomTopoTool, MB_CHK_SET_ERR, MB_SUCCESS, and moab::GeomTopoTool::setup_implicit_complement().

{

    ErrorCode 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;
}

ErrorCode moab::GeomQueryTool::measure_area	(	EntityHandle	surface,
		double &	result
	)

Calculate sum of area of triangles

Definition at line 1403 of file GeomQueryTool.cpp.

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.

{
    // 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;
}

ErrorCode moab::GeomQueryTool::measure_volume	(	EntityHandle	volume,
		double &	result
	)

Calculate the volume contained in a 'volume'

Definition at line 1340 of file GeomQueryTool.cpp.

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

{
    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;
}

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

Definition at line 381 of file GeomQueryTool.hpp.

References MBI.

    {
        return MBI;
    }

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 1081 of file GeomQueryTool.cpp.

References ErrorCode, geomTopoTool, moab::GeomTopoTool::get_bounding_coords(), and MB_CHK_SET_ERR.

Referenced by find_volume(), and point_in_volume().

{
    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;
}

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 920 of file GeomQueryTool.cpp.

References boundary_case(), counting, moab::debug, 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().

{
    // 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::cout << "direction==tangent" << std::endl;
                result = -1;
            }
            else
            {
                MB_SET_ERR( MB_FAILURE, "Error: unknown direction" );
            }
        }
    }

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

    return MB_SUCCESS;
}

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 1145 of file GeomQueryTool.cpp.

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

{
    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;
}

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 1564 of file GeomQueryTool.cpp.

References moab::CartVect::array(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::CartVect::length(), MB_CHK_SET_ERR, MB_SUCCESS, and MBI.

Referenced by point_in_volume_slow().

{
    ErrorCode 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;
}

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 739 of file GeomQueryTool.cpp.

References counting, moab::debug, 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.

{

    // 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;
        }
    }

    if( debug )
    {
        std::cout << "ray_fire:"
                  << " xyz=" << point[0] << " " << point[1] << " " << point[2] << " uvw=" << dir[0] << " " << dir[1]
                  << " " << dir[2] << " entity_handle=" << volume << std::endl;
    }

    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;
        if( debug )
        {
            std::cout << "          next_surf=0 dist=(undef)" << std::endl;
        }
        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;
        if( debug )
        {
            std::cout << "next surf hit = 0, dist = (undef)" << std::endl;
        }
        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] );
    }

    if( 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;
    }

    return MB_SUCCESS;
}

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 1631 of file GeomQueryTool.cpp.

References numericalPrecision.

{

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

    std::cout << "Set numerical precision = " << numericalPrecision << std::endl;
}

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 1617 of file GeomQueryTool.cpp.

References overlapThickness.

{

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

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 1107 of file GeomQueryTool.cpp.

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.

{
    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;
}