moab::AdaptiveKDTree Class Reference

Adaptive KD tree, for sorting and searching entities spatially. More...

#include <AdaptiveKDTree.hpp>

Inheritance diagram for moab::AdaptiveKDTree:

Collaboration diagram for moab::AdaptiveKDTree:

Classes
struct	Plane
	Split plane. More...

Public Types
enum	Axis { X = 0 , Y = 1 , Z = 2 }
	Enumeriate split plane directions. More...
enum	CandidatePlaneSet { SUBDIVISION = 0 , SUBDIVISION_SNAP , VERTEX_MEDIAN , VERTEX_SAMPLE }
	methods for selecting candidate split planes More...

Public Member Functions
	AdaptiveKDTree (Interface *iface)
	AdaptiveKDTree (Interface *iface, const Range &entities, EntityHandle *tree_root_set=NULL, FileOptions *opts=NULL)
	Constructor (build the tree on construction) Construct a tree object, and build the tree with entities input. See comments for build_tree() for detailed description of arguments. More...
	~AdaptiveKDTree ()
ErrorCode	parse_options (FileOptions &options)
	Parse options for tree creation. More...
virtual ErrorCode	build_tree (const Range &entities, EntityHandle *tree_root_set=NULL, FileOptions *options=NULL)
virtual ErrorCode	reset_tree ()
	Reset the tree, optionally checking we have the right root. More...
virtual ErrorCode	point_search (const double *point, EntityHandle &leaf_out, const double iter_tol=1.0e-10, const double inside_tol=1.0e-6, bool *multiple_leaves=NULL, EntityHandle *start_node=NULL, CartVect *params=NULL)
	Get leaf containing input position. More...
ErrorCode	point_search (const double *point, AdaptiveKDTreeIter &leaf_it, const double iter_tol=1.0e-10, const double inside_tol=1.0e-6, bool *multiple_leaves=NULL, EntityHandle *start_node=NULL)
	Get leaf containing input position. More...
virtual ErrorCode	distance_search (const double *point, const double distance, std::vector< EntityHandle > &leaves_out, const double iter_tol=1.0e-10, const double inside_tol=1.0e-6, std::vector< double > *dists_out=NULL, std::vector< CartVect > *params_out=NULL, EntityHandle *start_node=NULL)
	Find all leaves within a given distance from point If dists_out input non-NULL, also returns distances from each leaf; if point i is inside leaf, 0 is given as dists_out[i]. If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities in tree nodes and return containing entities in leaves_out. In those cases, if params_out is also non-NULL, will return parameters in those elements in that vector. More...
ErrorCode	get_info (EntityHandle root, double min[3], double max[3], unsigned int &dep)
	Return some basic information about the tree Stats are returned for tree starting from input node or tree root (root = 0) More...
ErrorCode	get_split_plane (EntityHandle node, Plane &plane)
	Get split plane for tree node. More...
ErrorCode	set_split_plane (EntityHandle node, const Plane &plane)
	Set split plane for tree node. More...
ErrorCode	get_tree_iterator (EntityHandle tree_root, AdaptiveKDTreeIter &result)
	Get iterator for tree. More...
ErrorCode	get_last_iterator (EntityHandle tree_root, AdaptiveKDTreeIter &result)
	Get iterator at right-most ('last') leaf. More...
ErrorCode	get_sub_tree_iterator (EntityHandle tree_root, const double box_min[3], const double box_max[3], AdaptiveKDTreeIter &result)
	Get iterator for tree or subtree. More...
ErrorCode	split_leaf (AdaptiveKDTreeIter &leaf, Plane plane)
	Split leaf of tree Updates iterator location to point to first new leaf node. More...
ErrorCode	split_leaf (AdaptiveKDTreeIter &leaf, Plane plane, EntityHandle &left_child, EntityHandle &right_child)
	Split leaf of tree Updates iterator location to point to first new leaf node. More...
ErrorCode	split_leaf (AdaptiveKDTreeIter &leaf, Plane plane, const Range &left_entities, const Range &right_entities)
	Split leaf of tree Updates iterator location to point to first new leaf node. More...
ErrorCode	split_leaf (AdaptiveKDTreeIter &leaf, Plane plane, const std::vector< EntityHandle > &left_entities, const std::vector< EntityHandle > &right_entities)
	Split leaf of tree Updates iterator location to point to first new leaf node. More...
ErrorCode	merge_leaf (AdaptiveKDTreeIter &iter)
	Merge the leaf pointed to by the current iterator with it's sibling. If the sibling is not a leaf, multiple merges may be done. More...
ErrorCode	closest_triangle (EntityHandle tree_root, const double from_coords[3], double closest_point_out[3], EntityHandle &triangle_out)
	Find triangle closest to input position. More...
ErrorCode	sphere_intersect_triangles (EntityHandle tree_root, const double center[3], double radius, std::vector< EntityHandle > &triangles)
ErrorCode	ray_intersect_triangles (EntityHandle tree_root, const double tolerance, const double ray_unit_dir[3], const double ray_base_pt[3], std::vector< EntityHandle > &triangles_out, std::vector< double > &distance_out, int result_count_limit=0, double distance_limit=-1.0)
ErrorCode	compute_depth (EntityHandle root, unsigned int &min_depth, unsigned int &max_depth)
virtual ErrorCode	print ()
	print various things about this tree More...
Public Member Functions inherited from moab::Tree
	Tree (Interface *iface)
	Constructor (bare) More...
virtual	~Tree ()
	Destructor. More...
ErrorCode	delete_tree_sets ()
	Delete the entity sets associated with the tree, starting with the root and traversing children. More...
virtual ErrorCode	get_bounding_box (BoundBox &box, EntityHandle *tree_node=NULL) const
	Get bounding box for tree below tree_node, or entire tree If no tree has been built yet, returns +/- DBL_MAX for all dimensions. Note for some tree types, boxes are not available for non-root nodes, and this function will return failure if non-root is passed in. More...
ErrorCode	find_all_trees (Range &results)
	Find all trees, by bounding box tag. More...
Interface *	moab ()
	Return the MOAB interface associated with this tree. More...
const Interface *	moab () const
	Return the MOAB interface associated with this tree. More...
double	get_max_depth ()
	Get max depth set on tree. More...
double	get_max_per_leaf ()
	Get max entities per leaf set on tree. More...
TreeStats &	tree_stats ()
	Get tree traversal stats object. More...
const TreeStats &	tree_stats () const
	Get tree traversal stats object. More...
ErrorCode	create_root (const double box_min[3], const double box_max[3], EntityHandle &root_handle)
	Create tree root and tag with bounding box. More...
ElemEvaluator *	get_eval ()
	get/set the ElemEvaluator More...
void	set_eval (ElemEvaluator *eval)
	get/set the ElemEvaluator More...

Private Member Functions
ErrorCode	init ()
ErrorCode	find_close_triangle (EntityHandle root, const double from_point[3], double pt[3], EntityHandle &triangle)
	find a triangle near the input point More...
ErrorCode	make_tag (Interface *iface, std::string name, TagType storage, DataType type, int count, void *default_val, Tag &tag_handle, std::vector< Tag > &created_tags)
ErrorCode	intersect_children_with_elems (const Range &elems, AdaptiveKDTree::Plane plane, double eps, CartVect box_min, CartVect box_max, Range &left_tris, Range &right_tris, Range &both_tris, double &metric_value)
ErrorCode	best_subdivision_snap_plane (int num_planes, const AdaptiveKDTreeIter &iter, Range &best_left, Range &best_right, Range &best_both, AdaptiveKDTree::Plane &best_plane, std::vector< double > &tmp_data, double eps)
ErrorCode	best_subdivision_plane (int num_planes, const AdaptiveKDTreeIter &iter, Range &best_left, Range &best_right, Range &best_both, AdaptiveKDTree::Plane &best_plane, double eps)
ErrorCode	best_vertex_median_plane (int num_planes, const AdaptiveKDTreeIter &iter, Range &best_left, Range &best_right, Range &best_both, AdaptiveKDTree::Plane &best_plane, std::vector< double > &coords, double eps)
ErrorCode	best_vertex_sample_plane (int num_planes, const AdaptiveKDTreeIter &iter, Range &best_left, Range &best_right, Range &best_both, AdaptiveKDTree::Plane &best_plane, std::vector< double > &coords, std::vector< EntityHandle > &indices, double eps)

Private Attributes
Tag	planeTag
Tag	axisTag
unsigned	splitsPerDir
CandidatePlaneSet	planeSet
bool	spherical
double	radius

Static Private Attributes
static const char *	treeName = "AKDTree"

Friends
class	AdaptiveKDTreeIter

Additional Inherited Members
Protected Member Functions inherited from moab::Tree
ErrorCode	parse_common_options (FileOptions &options)
	Parse options common to all trees. More...
Tag	get_box_tag (bool create_if_missing=true)
	Get the box tag, possibly constructing it first. More...
Protected Attributes inherited from moab::Tree
Interface *	mbImpl
BoundBox	boundBox
int	maxPerLeaf
int	maxDepth
int	treeDepth
double	minWidth
unsigned int	meshsetFlags
bool	cleanUp
EntityHandle	myRoot
Tag	boxTag
std::string	boxTagName
TreeStats	treeStats
ElemEvaluator *	myEval

Detailed Description

Adaptive KD tree, for sorting and searching entities spatially.

Examples

KDTree.cpp.

Definition at line 23 of file AdaptiveKDTree.hpp.

Member Enumeration Documentation

Axis

enum moab::AdaptiveKDTree::Axis

Enumeriate split plane directions.

Enumerator
X
Y
Z

Definition at line 148 of file AdaptiveKDTree.hpp.

    {
        X = 0,
        Y = 1,
        Z = 2
    };

CandidatePlaneSet

enum moab::AdaptiveKDTree::CandidatePlaneSet

methods for selecting candidate split planes

Enumerator
SUBDIVISION	Candidiate planes at evenly spaced intervals.
SUBDIVISION_SNAP	Like SUBDIVISION, except snap to closest vertex coordinate.
VERTEX_MEDIAN	Median vertex coodinate values.
VERTEX_SAMPLE	Random sampling of vertex coordinate values.

Definition at line 248 of file AdaptiveKDTree.hpp.

    {
        //! Candidiate planes at evenly spaced intervals
        SUBDIVISION = 0,
        //! Like SUBDIVISION, except snap to closest vertex coordinate
        SUBDIVISION_SNAP,  // = 1
                           //! Median vertex coodinate values
        VERTEX_MEDIAN,     // = 2
                           //! Random sampling of vertex coordinate values
        VERTEX_SAMPLE      // = 3
    };

Constructor & Destructor Documentation

AdaptiveKDTree() [1/2]

moab::AdaptiveKDTree::AdaptiveKDTree

(

Interface *

iface

)

Definition at line 36 of file AdaptiveKDTree.cpp.

    : Tree( iface ), planeTag( 0 ), axisTag( 0 ), splitsPerDir( 3 ), planeSet( SUBDIVISION_SNAP ), spherical( false ),
      radius( 1.0 )
{
    boxTagName = treeName;
 
    ErrorCode rval = init();
    if( MB_SUCCESS != rval ) throw rval;
}

References moab::Tree::boxTagName, ErrorCode, init(), MB_SUCCESS, and treeName.

AdaptiveKDTree() [2/2]

moab::AdaptiveKDTree::AdaptiveKDTree	(	Interface *	iface,
		const Range &	entities,
		EntityHandle *	tree_root_set = NULL,
		FileOptions *	opts = NULL
	)

Constructor (build the tree on construction) Construct a tree object, and build the tree with entities input. See comments for build_tree() for detailed description of arguments.

Parameters

iface	MOAB instance
entities	Entities to build tree around
tree_root	Root set for tree (see function description)
opts	Options for tree (see function description)

Definition at line 46 of file AdaptiveKDTree.cpp.

    : Tree( iface ), planeTag( 0 ), axisTag( 0 ), splitsPerDir( 3 ), planeSet( SUBDIVISION_SNAP ), spherical( false ),
      radius( 1.0 )
{
    boxTagName = treeName;
 
    ErrorCode rval;
    if( opts )
    {
        rval = parse_options( *opts );
        if( MB_SUCCESS != rval ) throw rval;
    }
 
    rval = init();
    if( MB_SUCCESS != rval ) throw rval;
 
    rval = build_tree( entities, tree_root_set, opts );
    if( MB_SUCCESS != rval ) throw rval;
}

References moab::Tree::boxTagName, build_tree(), ErrorCode, init(), MB_SUCCESS, parse_options(), and treeName.

~AdaptiveKDTree()

moab::AdaptiveKDTree::~AdaptiveKDTree

(

)

Definition at line 69 of file AdaptiveKDTree.cpp.

{
    if( !cleanUp ) return;
 
    if( myRoot )
    {
        reset_tree();
        myRoot = 0;
    }
}

References moab::Tree::cleanUp, moab::Tree::myRoot, and reset_tree().

Member Function Documentation

best_subdivision_plane()

ErrorCode moab::AdaptiveKDTree::best_subdivision_plane ( int  num_planes, const AdaptiveKDTreeIter &  iter, Range &  best_left, Range &  best_right, Range &  best_both, AdaptiveKDTree::Plane &  best_plane, double  eps  )

private

Definition at line 942 of file AdaptiveKDTree.cpp.

{
    double metric_val = std::numeric_limits< unsigned >::max();
 
    ErrorCode r;
    const CartVect box_min( iter.box_min() );
    const CartVect box_max( iter.box_max() );
    const CartVect diff( box_max - box_min );
 
    Range entities;
    r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
    if( MB_SUCCESS != r ) return r;
    const size_t p_count = entities.size();
 
    for( int axis = 0; axis < 3; ++axis )
    {
        int plane_count = num_planes;
        if( ( num_planes + 1 ) * eps >= diff[axis] ) plane_count = (int)( diff[axis] / eps ) - 1;
 
        for( int p = 1; p <= plane_count; ++p )
        {
            AdaptiveKDTree::Plane plane = { box_min[axis] + ( p / ( 1.0 + plane_count ) ) * diff[axis], axis };
            Range left, right, both;
            double val;
            r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
            if( MB_SUCCESS != r ) return r;
            const size_t sdiff = p_count - both.size();
            if( left.size() == sdiff || right.size() == sdiff ) continue;
 
            if( val >= metric_val ) continue;
 
            metric_val = val;
            best_plane = plane;
            best_left.swap( left );
            best_right.swap( right );
            best_both.swap( both );
        }
    }
 
    return MB_SUCCESS;
}

References moab::AdaptiveKDTreeIter::box_max(), box_max(), moab::AdaptiveKDTreeIter::box_min(), box_min(), ErrorCode, moab::Interface::get_entities_by_handle(), moab::AdaptiveKDTreeIter::handle(), intersect_children_with_elems(), MB_SUCCESS, moab::Tree::moab(), moab::Range::size(), moab::Range::swap(), and moab::AdaptiveKDTreeIter::tool().

Referenced by build_tree().

best_subdivision_snap_plane()

ErrorCode moab::AdaptiveKDTree::best_subdivision_snap_plane ( int  num_planes, const AdaptiveKDTreeIter &  iter, Range &  best_left, Range &  best_right, Range &  best_both, AdaptiveKDTree::Plane &  best_plane, std::vector< double > &  tmp_data, double  eps  )

private

Definition at line 990 of file AdaptiveKDTree.cpp.

{
    double metric_val = std::numeric_limits< unsigned >::max();
 
    ErrorCode r;
    // const CartVect tol(eps*diff);
 
    Range entities, vertices;
    r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
    if( MB_SUCCESS != r ) return r;
    const size_t p_count = entities.size();
    r                    = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
    if( MB_SUCCESS != r ) return r;
 
    unsigned int nverts = vertices.size();
    tmp_data.resize( 3 * nverts );
    r = iter.tool()->moab()->get_coords( vertices, &tmp_data[0], &tmp_data[nverts], &tmp_data[2 * nverts] );
    if( MB_SUCCESS != r ) return r;
 
    // calculate bounding box of vertices
    // decide based on the actual box the splitting plane
    // do not decide based on iterator box.
    // it could be too big
    // BoundBox box;
    // r = box.update(*moab(), vertices);
    CartVect box_min;
    CartVect box_max;
    for( int dir = 0; dir < 3; dir++ )
    {
        double amin = tmp_data[dir * nverts];
        double amax = amin;
        double* p   = &tmp_data[dir * nverts + 1];
        for( unsigned int i = 1; i < nverts; i++ )
        {
            if( *p < amin ) amin = *p;
            if( *p > amax ) amax = *p;
            p++;
        }
        box_min[dir] = amin;
        box_max[dir] = amax;
    }
    CartVect diff( box_max - box_min );
 
    for( int axis = 0; axis < 3; ++axis )
    {
        int plane_count = num_planes;
 
        // if num_planes results in width < eps, reset the plane count
        if( eps > 0.0 && std::isfinite(diff[axis]) && ( num_planes + 1 ) * eps >= diff[axis] )
          plane_count = (int)( diff[axis] / eps ) - 1;
 
        for( int p = 1; p <= plane_count; ++p )
        {
 
            // coord of this plane on axis
            double coord = box_min[axis] + ( p / ( 1.0 + plane_count ) ) * diff[axis];
 
            // find closest vertex coordinate to this plane position
            unsigned int istrt   = axis * nverts;
            double closest_coord = tmp_data[istrt];
            for( unsigned i = 1; i < nverts; ++i )
                if( fabs( coord - tmp_data[istrt + i] ) < fabs( coord - closest_coord ) )
                    closest_coord = tmp_data[istrt + i];
            if( closest_coord - box_min[axis] <= eps || box_max[axis] - closest_coord <= eps ) continue;
 
            // seprate elems into left/right/both, and compute separating metric
            AdaptiveKDTree::Plane plane = { closest_coord, axis };
            Range left, right, both;
            double val;
            r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
            if( MB_SUCCESS != r ) return r;
            const size_t d = p_count - both.size();
            if( left.size() == d || right.size() == d ) continue;
 
            if( val >= metric_val ) continue;
 
            metric_val = val;
            best_plane = plane;
            best_left.swap( left );
            best_right.swap( right );
            best_both.swap( both );
        }
    }
 
    return MB_SUCCESS;
}

References box_max(), box_min(), ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_coords(), moab::Interface::get_entities_by_handle(), moab::AdaptiveKDTreeIter::handle(), intersect_children_with_elems(), MB_SUCCESS, moab::Tree::moab(), moab::Range::size(), moab::Range::swap(), moab::AdaptiveKDTreeIter::tool(), and moab::Interface::UNION.

Referenced by build_tree().

best_vertex_median_plane()

ErrorCode moab::AdaptiveKDTree::best_vertex_median_plane ( int  num_planes, const AdaptiveKDTreeIter &  iter, Range &  best_left, Range &  best_right, Range &  best_both, AdaptiveKDTree::Plane &  best_plane, std::vector< double > &  coords, double  eps  )

private

Definition at line 1084 of file AdaptiveKDTree.cpp.

{
    double metric_val = std::numeric_limits< unsigned >::max();
 
    ErrorCode r;
    const CartVect box_min( iter.box_min() );
    const CartVect box_max( iter.box_max() );
 
    Range entities, vertices;
    r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
    if( MB_SUCCESS != r ) return r;
    const size_t p_count = entities.size();
    r                    = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
    if( MB_SUCCESS != r ) return r;
 
    coords.resize( vertices.size() );
    for( int axis = 0; axis < 3; ++axis )
    {
        if( box_max[axis] - box_min[axis] <= 2 * eps ) continue;
 
        double* ptrs[] = { 0, 0, 0 };
        ptrs[axis]     = &coords[0];
        r              = iter.tool()->moab()->get_coords( vertices, ptrs[0], ptrs[1], ptrs[2] );
        if( MB_SUCCESS != r ) return r;
 
        std::sort( coords.begin(), coords.end() );
        std::vector< double >::iterator citer;
        citer              = std::upper_bound( coords.begin(), coords.end(), box_min[axis] + eps );
        const size_t count = std::upper_bound( citer, coords.end(), box_max[axis] - eps ) - citer;
        size_t step;
        int np = num_planes;
        if( count < 2 * (size_t)num_planes )
        {
            step = 1;
            np   = count - 1;
        }
        else
        {
            step = count / ( num_planes + 1 );
        }
 
        for( int p = 1; p <= np; ++p )
        {
 
            citer += step;
            AdaptiveKDTree::Plane plane = { *citer, axis };
            Range left, right, both;
            double val;
            r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
            if( MB_SUCCESS != r ) return r;
            const size_t diff = p_count - both.size();
            if( left.size() == diff || right.size() == diff ) continue;
 
            if( val >= metric_val ) continue;
 
            metric_val = val;
            best_plane = plane;
            best_left.swap( left );
            best_right.swap( right );
            best_both.swap( both );
        }
    }
 
    return MB_SUCCESS;
}

References moab::AdaptiveKDTreeIter::box_max(), box_max(), moab::AdaptiveKDTreeIter::box_min(), box_min(), ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_coords(), moab::Interface::get_entities_by_handle(), moab::AdaptiveKDTreeIter::handle(), intersect_children_with_elems(), MB_SUCCESS, moab::Tree::moab(), moab::Range::size(), moab::Range::swap(), moab::AdaptiveKDTreeIter::tool(), and moab::Interface::UNION.

Referenced by build_tree().

best_vertex_sample_plane()

ErrorCode moab::AdaptiveKDTree::best_vertex_sample_plane ( int  num_planes, const AdaptiveKDTreeIter &  iter, Range &  best_left, Range &  best_right, Range &  best_both, AdaptiveKDTree::Plane &  best_plane, std::vector< double > &  coords, std::vector< EntityHandle > &  indices, double  eps  )

private

Definition at line 1157 of file AdaptiveKDTree.cpp.

{
    const size_t random_elem_threshold = 20 * num_planes;
    double metric_val                  = std::numeric_limits< unsigned >::max();
 
    ErrorCode r;
    const CartVect box_min( iter.box_min() );
    const CartVect box_max( iter.box_max() );
 
    Range entities, vertices;
    r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
    if( MB_SUCCESS != r ) return r;
 
    // We are selecting random vertex coordinates to use for candidate split
    // planes.  So if element list is large, begin by selecting random elements.
    const size_t p_count = entities.size();
    coords.resize( 3 * num_planes );
    if( p_count < random_elem_threshold )
    {
        r = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
        if( MB_SUCCESS != r ) return r;
    }
    else
    {
        indices.resize( random_elem_threshold );
        const int num_rand = p_count / RAND_MAX + 1;
        for( size_t j = 0; j < random_elem_threshold; ++j )
        {
            size_t rnd = rand();
            for( int i = num_rand; i > 1; --i )
                rnd *= rand();
            rnd %= p_count;
            indices[j] = entities[rnd];
        }
        r = iter.tool()->moab()->get_adjacencies( &indices[0], random_elem_threshold, 0, false, vertices,
                                                  Interface::UNION );
        if( MB_SUCCESS != r ) return r;
    }
 
    coords.resize( vertices.size() );
    for( int axis = 0; axis < 3; ++axis )
    {
        if( box_max[axis] - box_min[axis] <= 2 * eps ) continue;
 
        double* ptrs[] = { 0, 0, 0 };
        ptrs[axis]     = &coords[0];
        r              = iter.tool()->moab()->get_coords( vertices, ptrs[0], ptrs[1], ptrs[2] );
        if( MB_SUCCESS != r ) return r;
 
        size_t num_valid_coords = 0;
        for( size_t i = 0; i < coords.size(); ++i )
            if( coords[i] > box_min[axis] + eps && coords[i] < box_max[axis] - eps ) ++num_valid_coords;
 
        if( 2 * (size_t)num_planes > num_valid_coords )
        {
            indices.clear();
            for( size_t i = 0; i < coords.size(); ++i )
                if( coords[i] > box_min[axis] + eps && coords[i] < box_max[axis] - eps ) indices.push_back( i );
        }
        else
        {
            indices.resize( num_planes );
            // make sure random indices are sufficient to cover entire range
            const int num_rand = coords.size() / RAND_MAX + 1;
            for( int j = 0; j < num_planes; ++j )
            {
                size_t rnd;
                do
                {
                    rnd = rand();
                    for( int i = num_rand; i > 1; --i )
                        rnd *= rand();
                    rnd %= coords.size();
                } while( coords[rnd] <= box_min[axis] + eps || coords[rnd] >= box_max[axis] - eps );
                indices[j] = rnd;
            }
        }
 
        for( unsigned p = 0; p < indices.size(); ++p )
        {
 
            AdaptiveKDTree::Plane plane = { coords[indices[p]], axis };
            Range left, right, both;
            double val;
            r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
            if( MB_SUCCESS != r ) return r;
            const size_t diff = p_count - both.size();
            if( left.size() == diff || right.size() == diff ) continue;
 
            if( val >= metric_val ) continue;
 
            metric_val = val;
            best_plane = plane;
            best_left.swap( left );
            best_right.swap( right );
            best_both.swap( both );
        }
    }
 
    return MB_SUCCESS;
}

References moab::AdaptiveKDTreeIter::box_max(), box_max(), moab::AdaptiveKDTreeIter::box_min(), box_min(), ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_coords(), moab::Interface::get_entities_by_handle(), moab::AdaptiveKDTreeIter::handle(), intersect_children_with_elems(), MB_SUCCESS, moab::Tree::moab(), moab::Range::size(), moab::Range::swap(), moab::AdaptiveKDTreeIter::tool(), and moab::Interface::UNION.

Referenced by build_tree().

build_tree()

ErrorCode moab::AdaptiveKDTree::build_tree ( const Range &  entities, EntityHandle *  tree_root_set = NULL, FileOptions *  options = NULL  )

virtual

Build the tree Build a tree with the entities input. If a non-NULL tree_root_set pointer is input, use the pointed-to set as the root of this tree (*tree_root_set!=0) otherwise construct a new root set and pass its handle back in *tree_root_set. Options vary by tree type; see Tree.hpp for common options; options specific to AdaptiveKDTree: SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3 PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below) for possible values; default = 1 (SUBDIVISION_SNAP)

Parameters

entities	Entities with which to build the tree
tree_root	Root set for tree (see function description)
opts	Options for tree (see function description)

Returns

Error is returned only on build failure

Implements moab::Tree.

Definition at line 80 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval;
    CpuTimer cp;
 
    if( options )
    {
        rval = parse_options( *options );
        if( MB_SUCCESS != rval ) return rval;
 
        if( !options->all_seen() ) return MB_FAILURE;
    }
 
    // calculate bounding box of elements
    BoundBox box;
    rval = box.update( *moab(), entities, spherical, radius );
    if( MB_SUCCESS != rval ) return rval;
 
    // create tree root
    EntityHandle tmp_root;
    if( !tree_root_set ) tree_root_set = &tmp_root;
    rval = create_root( box.bMin.array(), box.bMax.array(), *tree_root_set );
    if( MB_SUCCESS != rval ) return rval;
    rval = moab()->add_entities( *tree_root_set, entities );
    if( MB_SUCCESS != rval ) return rval;
 
    AdaptiveKDTreeIter iter;
    iter.initialize( this, *tree_root_set, box.bMin.array(), box.bMax.array(), AdaptiveKDTreeIter::LEFT );
 
    std::vector< double > tmp_data;
    std::vector< EntityHandle > tmp_data2;
    for( ;; )
    {
 
        int pcount;
        rval = moab()->get_number_entities_by_handle( iter.handle(), pcount );
        if( MB_SUCCESS != rval ) break;
 
        const size_t p_count = pcount;
        Range best_left, best_right, best_both;
        Plane best_plane = { HUGE_VAL, -1 };
        if( (int)p_count > maxPerLeaf && (int)iter.depth() < maxDepth )
        {
            switch( planeSet )
            {
                case AdaptiveKDTree::SUBDIVISION:
                    rval = best_subdivision_plane( splitsPerDir, iter, best_left, best_right, best_both, best_plane,
                                                   minWidth );
                    break;
                case AdaptiveKDTree::SUBDIVISION_SNAP:
                    rval = best_subdivision_snap_plane( splitsPerDir, iter, best_left, best_right, best_both,
                                                        best_plane, tmp_data, minWidth );
                    break;
                case AdaptiveKDTree::VERTEX_MEDIAN:
                    rval = best_vertex_median_plane( splitsPerDir, iter, best_left, best_right, best_both, best_plane,
                                                     tmp_data, minWidth );
                    break;
                case AdaptiveKDTree::VERTEX_SAMPLE:
                    rval = best_vertex_sample_plane( splitsPerDir, iter, best_left, best_right, best_both, best_plane,
                                                     tmp_data, tmp_data2, minWidth );
                    break;
                default:
                    rval = MB_FAILURE;
            }
 
            if( MB_SUCCESS != rval ) return rval;
        }
 
        if( best_plane.norm >= 0 )
        {
            best_left.merge( best_both );
            best_right.merge( best_both );
            rval = split_leaf( iter, best_plane, best_left, best_right );
            if( MB_SUCCESS != rval ) return rval;
        }
        else
        {
            rval = iter.step();
            if( MB_ENTITY_NOT_FOUND == rval )
            {
                rval               = treeStats.compute_stats( mbImpl, myRoot );
                treeStats.initTime = cp.time_elapsed();
                return rval;  // at end
            }
            else if( MB_SUCCESS != rval )
                break;
        }
    }
 
    reset_tree();
 
    treeStats.reset();
 
    return rval;
}

References moab::Interface::add_entities(), moab::FileOptions::all_seen(), moab::CartVect::array(), best_subdivision_plane(), best_subdivision_snap_plane(), best_vertex_median_plane(), best_vertex_sample_plane(), moab::BoundBox::bMax, moab::BoundBox::bMin, moab::TreeStats::compute_stats(), moab::Tree::create_root(), moab::AdaptiveKDTreeIter::depth(), ErrorCode, moab::Interface::get_number_entities_by_handle(), moab::AdaptiveKDTreeIter::handle(), moab::AdaptiveKDTreeIter::initialize(), moab::TreeStats::initTime, moab::AdaptiveKDTreeIter::LEFT, moab::Tree::maxDepth, moab::Tree::maxPerLeaf, MB_ENTITY_NOT_FOUND, MB_SUCCESS, moab::Tree::mbImpl, moab::Range::merge(), moab::Tree::minWidth, moab::Tree::moab(), moab::Tree::myRoot, moab::AdaptiveKDTree::Plane::norm, parse_options(), planeSet, radius, moab::TreeStats::reset(), reset_tree(), spherical, split_leaf(), splitsPerDir, moab::AdaptiveKDTreeIter::step(), SUBDIVISION, SUBDIVISION_SNAP, moab::CpuTimer::time_elapsed(), moab::Tree::treeStats, moab::BoundBox::update(), VERTEX_MEDIAN, and VERTEX_SAMPLE.

Referenced by AdaptiveKDTree(), moab::TempestRemapper::ConstructCoveringSet(), moab::Coupler::initialize_tree(), moab::Intx2Mesh::intersect_meshes(), moab::Intx2Mesh::intersect_meshes_kdtree(), main(), moab::MergeMesh::merge_all(), merge_duplicate_vertices(), moab::MergeMesh::merge_entities(), and moab::ReadNCDF::update().

closest_triangle()

ErrorCode moab::AdaptiveKDTree::closest_triangle	(	EntityHandle	tree_root,
		const double	from_coords[3],
		double	closest_point_out[3],
		EntityHandle &	triangle_out
	)

Find triangle closest to input position.

Parameters

from_coords	The input position to test against
closest_point_out	The closest point on the set of triangles in the tree
triangle_out	The triangle closest to the input position

Find the triangles in a set that are closer to the input position than any triangles in the 'closest_tris' list.

closest_tris is assumed to contain a list of triangles for which the first is the closest known triangle to the input position and the first entry in 'closest_pts' is the closest location on that triangle. Any other values in the lists must be other triangles for which the closest point is within the input tolerance of the closest closest point. This function will update the lists as appropriate if any closer triangles or triangles within the tolerance of the current closest location are found. The first entry is maintained as the closest of the list of triangles.

Definition at line 1764 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval;
    double shortest_dist_sqr = HUGE_VAL;
    std::vector< EntityHandle > leaves;
    const CartVect from( from_coords );
    CartVect closest_pt;
 
    // Find the leaf containing the input point
    // This search does not take into account any bounding box for the
    // tree, so it always returns one leaf.
    assert( tree_root );
    rval = find_close_triangle( tree_root, from_coords, closest_pt.array(), triangle_out );
    if( MB_SUCCESS != rval ) return rval;
 
    // Find any other leaves for which the bounding box is within
    // the same distance from the input point as the current closest
    // point is.
    CartVect diff = closest_pt - from;
    rval = distance_search( from_coords, sqrt( diff % diff ), leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root );
    if( MB_SUCCESS != rval ) return rval;
 
    // Check any close leaves to see if they contain triangles that
    // are as close to or closer than the current closest triangle(s).
    for( unsigned i = 0; i < leaves.size(); ++i )
    {
        rval = closest_to_triangles( moab(), leaves[i], from, shortest_dist_sqr, closest_pt, triangle_out );
        if( MB_SUCCESS != rval ) return rval;
    }
 
    // pass back resulting position
    closest_pt.get( closest_point_out );
    return MB_SUCCESS;
}

References moab::CartVect::array(), moab::closest_to_triangles(), distance_search(), ErrorCode, find_close_triangle(), moab::CartVect::get(), MB_SUCCESS, and moab::Tree::moab().

compute_depth()

ErrorCode moab::AdaptiveKDTree::compute_depth	(	EntityHandle	root,
		unsigned int &	min_depth,
		unsigned int &	max_depth
	)

Definition at line 2032 of file AdaptiveKDTree.cpp.

{
    AdaptiveKDTreeIter iter;
    get_tree_iterator( root, iter );
    iter.step_to_first_leaf( AdaptiveKDTreeIter::LEFT );
    min_depth = max_depth = iter.depth();
 
    int num_of_elements = 0, max, min;
    moab()->get_number_entities_by_handle( iter.handle(), num_of_elements );
    max = min = num_of_elements;
    while( MB_SUCCESS == iter.step() )
    {
        int temp = 0;
        moab()->get_number_entities_by_handle( iter.handle(), temp );
        max = std::max( max, temp );
        min = std::min( min, temp );
        if( iter.depth() > max_depth )
            max_depth = iter.depth();
        else if( iter.depth() < min_depth )
            min_depth = iter.depth();
    }
    return MB_SUCCESS;
}

References moab::AdaptiveKDTreeIter::depth(), moab::Interface::get_number_entities_by_handle(), get_tree_iterator(), moab::AdaptiveKDTreeIter::handle(), moab::AdaptiveKDTreeIter::LEFT, MB_SUCCESS, moab::Tree::moab(), moab::AdaptiveKDTreeIter::step(), and moab::AdaptiveKDTreeIter::step_to_first_leaf().

Referenced by get_info().

distance_search()

ErrorCode moab::AdaptiveKDTree::distance_search ( const double *  point, const double  distance, std::vector< EntityHandle > &  leaves_out, const double  iter_tol = 1.0e-10, const double  inside_tol = 1.0e-6, std::vector< double > *  dists_out = NULL, std::vector< CartVect > *  params_out = NULL, EntityHandle *  start_node = NULL  )

virtual

Find all leaves within a given distance from point If dists_out input non-NULL, also returns distances from each leaf; if point i is inside leaf, 0 is given as dists_out[i]. If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities in tree nodes and return containing entities in leaves_out. In those cases, if params_out is also non-NULL, will return parameters in those elements in that vector.

Parameters

point	Point to be located in tree
distance	Distance within which to query
leaves_out	Leaves within distance or containing point
iter_tol	Tolerance for convergence of point search
inside_tol	Tolerance for inside element calculation
dists_out	If non-NULL, will contain distsances to leaves
params_out	If non-NULL, will contain parameters of the point in the ents in leaves_out
start_node	Start from this tree node (non-NULL) instead of tree root (NULL)

when we pop it from the list.)

when we pop it from the list.)

Implements moab::Tree.

Definition at line 1389 of file AdaptiveKDTree.cpp.

{
    treeStats.numTraversals++;
    const double dist_sqr = distance * distance;
    const CartVect from( from_point );
    std::vector< NodeDistance > list,
        result_list_nodes;  // list of subtrees to traverse, and results
                            // pre-allocate space for default max tree depth
    list.reserve( maxDepth );
 
    // misc temporary values
    Plane plane;
    NodeDistance node;
    ErrorCode rval;
    std::vector< EntityHandle > children;
 
    // Get distance from input position to bounding box of tree
    // (zero if inside box)
    BoundBox box;
    rval = get_bounding_box( box );
    if( MB_SUCCESS == rval && !box.contains_point( from_point, iter_tol ) )
    {
        treeStats.nodesVisited++;
        return MB_SUCCESS;
    }
 
    // if bounding box is not available (e.g. not starting from true root)
    // just start with zero.  Less efficient, but will work.
    node.dist = CartVect( 0.0 );
    if( MB_SUCCESS == rval )
    {
        for( int i = 0; i < 3; ++i )
        {
            if( from_point[i] < box.bMin[i] )
                node.dist[i] = box.bMin[i] - from_point[i];
            else if( from_point[i] > box.bMax[i] )
                node.dist[i] = from_point[i] - box.bMax[i];
        }
        if( node.dist % node.dist > dist_sqr )
        {
            treeStats.nodesVisited++;
            return MB_SUCCESS;
        }
    }
 
    // begin with root in list
    node.handle = ( tree_root ? *tree_root : myRoot );
    list.push_back( node );
 
    while( !list.empty() )
    {
 
        node = list.back();
        list.pop_back();
        treeStats.nodesVisited++;
 
        // If leaf node, test contained triangles
        children.clear();
        rval = moab()->get_child_meshsets( node.handle, children );
        if( children.empty() )
        {
            treeStats.leavesVisited++;
            if( myEval && result_params )
            {
                EntityHandle ent;
                CartVect params;
                rval = myEval->find_containing_entity( node.handle, from_point, iter_tol, inside_tol, ent,
                                                       params.array(), &treeStats.traversalLeafObjectTests );
                if( MB_SUCCESS != rval )
                    return rval;
                else if( ent )
                {
                    result_list.push_back( ent );
                    result_params->push_back( params );
                    if( result_dists ) result_dists->push_back( 0.0 );
                }
            }
            else
            {
                result_list_nodes.push_back( node );
                continue;
            }
        }
 
        // If not leaf node, add children to working list
        rval = get_split_plane( node.handle, plane );
        if( MB_SUCCESS != rval ) return rval;
 
        const double d = from[plane.norm] - plane.coord;
 
        // right of plane?
        if( d > 0 )
        {
            node.handle = children[1];
            list.push_back( node );
            // if the split plane is close to the input point, add
            // the left child also (we'll check the exact distance
            /// when we pop it from the list.)
            if( d <= distance )
            {
                node.dist[plane.norm] = d;
                if( node.dist % node.dist <= dist_sqr )
                {
                    node.handle = children[0];
                    list.push_back( node );
                }
            }
        }
        // left of plane
        else
        {
            node.handle = children[0];
            list.push_back( node );
            // if the split plane is close to the input point, add
            // the right child also (we'll check the exact distance
            /// when we pop it from the list.)
            if( -d <= distance )
            {
                node.dist[plane.norm] = -d;
                if( node.dist % node.dist <= dist_sqr )
                {
                    node.handle = children[1];
                    list.push_back( node );
                }
            }
        }
    }
 
    if( myEval && result_params ) return MB_SUCCESS;
 
    // separate loops to avoid if test inside loop
 
    result_list.reserve( result_list_nodes.size() );
    for( std::vector< NodeDistance >::iterator vit = result_list_nodes.begin(); vit != result_list_nodes.end(); ++vit )
        result_list.push_back( ( *vit ).handle );
 
    if( result_dists && distance > 0.0 )
    {
        result_dists->reserve( result_list_nodes.size() );
        for( std::vector< NodeDistance >::iterator vit = result_list_nodes.begin(); vit != result_list_nodes.end();
             ++vit )
            result_dists->push_back( ( *vit ).dist.length() );
    }
 
    return MB_SUCCESS;
}

References moab::BoundBox::bMax, moab::BoundBox::bMin, moab::BoundBox::contains_point(), moab::AdaptiveKDTree::Plane::coord, moab::NodeDistance::dist, ErrorCode, moab::ElemEvaluator::find_containing_entity(), moab::Tree::get_bounding_box(), moab::Interface::get_child_meshsets(), get_split_plane(), moab::NodeDistance::handle, moab::TreeStats::leavesVisited, moab::Tree::maxDepth, MB_SUCCESS, moab::Tree::moab(), moab::Tree::myEval, moab::Tree::myRoot, moab::TreeStats::nodesVisited, moab::AdaptiveKDTree::Plane::norm, moab::TreeStats::numTraversals, moab::params, moab::TreeStats::traversalLeafObjectTests, and moab::Tree::treeStats.

Referenced by closest_triangle(), moab::MergeMesh::find_merged_to(), moab::Intx2Mesh::intersect_meshes(), moab::Intx2Mesh::intersect_meshes_kdtree(), merge_duplicate_vertices(), moab::Coupler::nat_param(), sphere_intersect_triangles(), and moab::ReadNCDF::update().

find_close_triangle()

ErrorCode moab::AdaptiveKDTree::find_close_triangle ( EntityHandle  root, const double  from_point[3], double  pt[3], EntityHandle &  triangle  )

private

find a triangle near the input point

Definition at line 1594 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval;
    Range tris;
    Plane split;
    std::vector< EntityHandle > stack;
    std::vector< EntityHandle > children( 2 );
    stack.reserve( 30 );
    assert( root );
    stack.push_back( root );
 
    while( !stack.empty() )
    {
        EntityHandle node = stack.back();
        stack.pop_back();
 
        for( ;; )
        {  // loop until we find a leaf
 
            children.clear();
            rval = moab()->get_child_meshsets( node, children );
            if( MB_SUCCESS != rval ) return rval;
 
            // loop termination criterion
            if( children.empty() ) break;
 
            // if not a leaf, get split plane
            rval = get_split_plane( node, split );
            if( MB_SUCCESS != rval ) return rval;
 
            // continue down the side that contains the point,
            // and push the other side onto the stack in case
            // we need to check it later.
            int rs = split.right_side( from );
            node   = children[rs];
            stack.push_back( children[1 - rs] );
        }
 
        // We should now be at a leaf.
        // If it has some triangles, we're done.
        // If not, continue searching for another leaf.
        tris.clear();
        rval = moab()->get_entities_by_type( node, MBTRI, tris );
        if( !tris.empty() )
        {
            double dist_sqr = HUGE_VAL;
            CartVect point( pt );
            rval = closest_to_triangles( moab(), tris, CartVect( from ), dist_sqr, point, triangle );
            point.get( pt );
            return rval;
        }
    }
 
    // If we got here, then we traversed the entire tree
    // and all the leaves were empty.
    return MB_ENTITY_NOT_FOUND;
}

References moab::Range::clear(), moab::closest_to_triangles(), moab::Range::empty(), ErrorCode, moab::CartVect::get(), moab::Interface::get_child_meshsets(), moab::Interface::get_entities_by_type(), get_split_plane(), MB_ENTITY_NOT_FOUND, MB_SUCCESS, MBTRI, moab::Tree::moab(), and split().

Referenced by closest_triangle().

get_info()

ErrorCode moab::AdaptiveKDTree::get_info ( EntityHandle  root, double  min[3], double  max[3], unsigned int &  max_dep  )

virtual

Return some basic information about the tree Stats are returned for tree starting from input node or tree root (root = 0)

Parameters

root	If non-0, give stats below and including root
min	Minimum corner of bounding box
max	Maximum corner of bounding box
max_dep	Maximum depth of tree below root

Reimplemented from moab::Tree.

Definition at line 2056 of file AdaptiveKDTree.cpp.

{
    BoundBox box;
    ErrorCode result = get_bounding_box( box, &root );
    if( MB_SUCCESS != result ) return result;
    box.bMin.get( bmin );
    box.bMax.get( bmax );
 
    unsigned min_depth;
    return compute_depth( root, min_depth, dep );
}

References moab::BoundBox::bMax, moab::BoundBox::bMin, compute_depth(), ErrorCode, moab::CartVect::get(), moab::Tree::get_bounding_box(), and MB_SUCCESS.

Referenced by moab::Coupler::initialize_tree().

get_last_iterator()

ErrorCode moab::AdaptiveKDTree::get_last_iterator	(	EntityHandle	tree_root,
		AdaptiveKDTreeIter &	result
	)

Get iterator at right-most ('last') leaf.

Definition at line 327 of file AdaptiveKDTree.cpp.

{
    double box[6];
    ErrorCode rval = moab()->tag_get_data( boxTag, &root, 1, box );
    if( MB_SUCCESS != rval ) return rval;
 
    return iter.initialize( this, root, box, box + 3, AdaptiveKDTreeIter::RIGHT );
}

References moab::Tree::boxTag, ErrorCode, moab::AdaptiveKDTreeIter::initialize(), MB_SUCCESS, moab::Tree::moab(), moab::AdaptiveKDTreeIter::RIGHT, and moab::Interface::tag_get_data().

get_split_plane()

ErrorCode moab::AdaptiveKDTree::get_split_plane	(	EntityHandle	node,
		Plane &	plane
	)

Get split plane for tree node.

Definition at line 285 of file AdaptiveKDTree.cpp.

{
#ifndef MB_AD_KD_TREE_USE_SINGLE_TAG
    ErrorCode r1, r2;
    r1 = moab()->tag_get_data( planeTag, &entity, 1, &plane.coord );
    r2 = moab()->tag_get_data( axisTag, &entity, 1, &plane.norm );
    return MB_SUCCESS == r1 ? r2 : r1;
#elif defined( MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG )
    double values[2];
    ErrorCode rval = moab()->tag_get_data( planeTag, &entity, 1, values );
    plane.coord    = values[0];
    plane.norm     = (int)values[1];
    return rval;
#else
    return moab()->tag_get_data( planeTag, &entity, 1, &plane );
#endif
}

References axisTag, moab::AdaptiveKDTree::Plane::coord, ErrorCode, MB_SUCCESS, moab::Tree::moab(), moab::AdaptiveKDTree::Plane::norm, planeTag, and moab::Interface::tag_get_data().

Referenced by distance_search(), find_close_triangle(), moab::AdaptiveKDTreeIter::get_neighbors(), moab::AdaptiveKDTreeIter::get_parent_split_plane(), merge_leaf(), point_search(), ray_intersect_triangles(), moab::AdaptiveKDTreeIter::sibling_side(), moab::AdaptiveKDTreeIter::step(), and moab::AdaptiveKDTreeIter::step_to_first_leaf().

get_sub_tree_iterator()

ErrorCode moab::AdaptiveKDTree::get_sub_tree_iterator	(	EntityHandle	tree_root,
		const double	box_min[3],
		const double	box_max[3],
		AdaptiveKDTreeIter &	result
	)

Get iterator for tree or subtree.

Definition at line 336 of file AdaptiveKDTree.cpp.

{
    return result.initialize( this, root, min, max, AdaptiveKDTreeIter::LEFT );
}

References moab::AdaptiveKDTreeIter::initialize(), and moab::AdaptiveKDTreeIter::LEFT.

Referenced by get_tree_iterator().

get_tree_iterator()

ErrorCode moab::AdaptiveKDTree::get_tree_iterator	(	EntityHandle	tree_root,
		AdaptiveKDTreeIter &	result
	)

Get iterator for tree.

Definition at line 318 of file AdaptiveKDTree.cpp.

{
    double box[6];
    ErrorCode rval = moab()->tag_get_data( boxTag, &root, 1, box );
    if( MB_SUCCESS != rval ) return rval;
 
    return get_sub_tree_iterator( root, box, box + 3, iter );
}

References moab::Tree::boxTag, ErrorCode, get_sub_tree_iterator(), MB_SUCCESS, moab::Tree::moab(), and moab::Interface::tag_get_data().

Referenced by compute_depth(), moab::MergeMesh::find_merged_to(), moab::Coupler::nat_param(), print(), and moab::ReadNCDF::update().

init()

ErrorCode moab::AdaptiveKDTree::init ( )

private

Definition at line 240 of file AdaptiveKDTree.cpp.

{
    std::vector< Tag > ctl;
 
#ifndef MB_AD_KD_TREE_USE_SINGLE_TAG
    // create two tags, one for axis direction and one for axis coordinate
    std::string n1( treeName ), n2( treeName );
    n1 += "_coord";
    n2 += "_norm";
    ErrorCode rval = make_tag( moab(), n1, MB_TAG_DENSE, MB_TYPE_DOUBLE, 1, 0, planeTag, ctl );
    if( MB_SUCCESS != rval ) return rval;
    rval = make_tag( moab(), n2, MB_TAG_DENSE, MB_TYPE_INT, 1, 0, axisTag, ctl );
    if( MB_SUCCESS != rval ) return rval;
 
#elif defined( MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG )
    // create tag to hold two doubles, one for location and one for axis
    std::string double_tag_name = std::string( treeName ) + std::string( "_coord_norm" );
    ErrorCode rval = make_tag( moab(), double_tag_name, MB_TAG_DENSE, MB_TYPE_DOUBLE, 2, 0, planeTag, ctl );
    if( MB_SUCCESS != rval ) return rval;
#else
    // create opaque tag to hold struct Plane
    ErrorCode rval = make_tag( moab(), tagname, MB_TAG_DENSE, MB_TYPE_OPAQUE, sizeof( Plane ), 0, planeTag, ctl );
    if( MB_SUCCESS != rval ) return rval;
 
#ifdef MOAB_HAVE_HDF5
    // create a mesh tag holding the HDF5 type for a struct Plane
    Tag type_tag;
    std::string type_tag_name = "__hdf5_tag_type_";
    type_tag_name += boxTagName;
    rval = make_tag( moab(), type_tag_name, MB_TAG_MESH, MB_TYPE_OPAQUE, sizeof( hid_t ), 0, type_tag, ctl );
    if( MB_SUCCESS != rval ) return rval;
    // create HDF5 type object describing struct Plane
    Plane p;
    hid_t handle = H5Tcreate( H5T_COMPOUND, sizeof( Plane ) );
    H5Tinsert( handle, "coord", &( p.coord ) - &p, H5T_NATIVE_DOUBLE );
    H5Tinsert( handle, "norm", &( p.axis ) - &p, H5T_NATIVE_INT );
    EntityHandle root = 0;
    rval              = mbImpl->tag_set_data( type_tag, &root, 1, &handle );
    if( MB_SUCCESS != rval ) return rval;
#endif
#endif
 
    return rval;
}

References axisTag, moab::Tree::boxTagName, moab::AdaptiveKDTree::Plane::coord, ErrorCode, make_tag(), MB_SUCCESS, MB_TAG_DENSE, MB_TAG_MESH, MB_TYPE_DOUBLE, MB_TYPE_OPAQUE, moab::Tree::mbImpl, moab::Tree::moab(), planeTag, moab::Interface::tag_set_data(), and treeName.

Referenced by AdaptiveKDTree().

intersect_children_with_elems()

ErrorCode moab::AdaptiveKDTree::intersect_children_with_elems ( const Range &  elems, AdaptiveKDTree::Plane  plane, double  eps, CartVect  box_min, CartVect  box_max, Range &  left_tris, Range &  right_tris, Range &  both_tris, double &  metric_value  )

private

Definition at line 777 of file AdaptiveKDTree.cpp.

{
    left_tris.clear();
    right_tris.clear();
    both_tris.clear();
    CartVect coords[16];
 
    // get extents of boxes for left and right sides
    BoundBox left_box( box_min, box_max ), right_box( box_min, box_max );
    right_box.bMin             = box_min;
    left_box.bMax              = box_max;
    right_box.bMin[plane.norm] = left_box.bMax[plane.norm] = plane.coord;
    const CartVect left_cen                                = 0.5 * ( left_box.bMax + box_min );
    const CartVect left_dim                                = 0.5 * ( left_box.bMax - box_min );
    const CartVect right_cen                               = 0.5 * ( box_max + right_box.bMin );
    const CartVect right_dim                               = 0.5 * ( box_max - right_box.bMin );
    const CartVect dim                                     = box_max - box_min;
    const double max_tol                                   = std::max( dim[0], std::max( dim[1], dim[2] ) ) / 10;
 
    // test each entity
    ErrorCode rval;
    int count, count2;
    const EntityHandle *conn, *conn2;
 
    const Range::const_iterator elem_begin = elems.lower_bound( MBEDGE );
    const Range::const_iterator poly_begin = elems.lower_bound( MBPOLYHEDRON, elem_begin );
    const Range::const_iterator set_begin  = elems.lower_bound( MBENTITYSET, poly_begin );
    Range::iterator left_ins               = left_tris.begin();
    Range::iterator right_ins              = right_tris.begin();
    Range::iterator both_ins               = both_tris.begin();
    Range::const_iterator i;
 
    // vertices
    for( i = elems.begin(); i != elem_begin; ++i )
    {
        tree_stats().constructLeafObjectTests++;
        rval = moab()->get_coords( &*i, 1, coords[0].array() );
        if( MB_SUCCESS != rval ) return rval;
 
        bool lo = false, ro = false;
        if( coords[0][plane.norm] <= plane.coord ) lo = true;
        if( coords[0][plane.norm] >= plane.coord ) ro = true;
 
        if( lo && ro )
            both_ins = both_tris.insert( both_ins, *i, *i );
        else if( lo )
            left_ins = left_tris.insert( left_ins, *i, *i );
        else  // if (ro)
            right_ins = right_tris.insert( right_ins, *i, *i );
    }
 
    // non-polyhedron elements
    std::vector< EntityHandle > dum_vector;
    for( i = elem_begin; i != poly_begin; ++i )
    {
        tree_stats().constructLeafObjectTests++;
        rval = moab()->get_connectivity( *i, conn, count, true, &dum_vector );
        if( MB_SUCCESS != rval ) return rval;
        if( count > (int)( sizeof( coords ) / sizeof( coords[0] ) ) ) return MB_FAILURE;
        rval = moab()->get_coords( &conn[0], count, coords[0].array() );
        if( MB_SUCCESS != rval ) return rval;
 
        bool lo = false, ro = false;
        for( int j = 0; j < count; ++j )
        {
            if( coords[j][plane.norm] <= plane.coord ) lo = true;
            if( coords[j][plane.norm] >= plane.coord ) ro = true;
        }
 
        // Triangle must be in at least one leaf.  If test against plane
        // identified that leaf, then we're done.  If triangle is on both
        // sides of plane, do more precise test to ensure that it is really
        // in both.
        //        BoundBox box;
        //        box.update(*moab(), *i);
        if( lo && ro )
        {
            double tol = eps;
            lo = ro = false;
            while( !lo && !ro && tol <= max_tol )
            {
                tree_stats().boxElemTests += 2;
                lo = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE( *i ), left_cen, left_dim + CartVect( tol ),
                                                 count );
                ro = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE( *i ), right_cen, right_dim + CartVect( tol ),
                                                 count );
 
                tol *= 10.0;
            }
        }
        if( lo && ro )
            both_ins = both_tris.insert( both_ins, *i, *i );
        else if( lo )
            left_ins = left_tris.insert( left_ins, *i, *i );
        else if( ro )
            right_ins = right_tris.insert( right_ins, *i, *i );
    }
 
    // polyhedra
    for( i = poly_begin; i != set_begin; ++i )
    {
        tree_stats().constructLeafObjectTests++;
        rval = moab()->get_connectivity( *i, conn, count, true );
        if( MB_SUCCESS != rval ) return rval;
 
        // just check the bounding box of the polyhedron
        bool lo = false, ro = false;
        for( int j = 0; j < count; ++j )
        {
            rval = moab()->get_connectivity( conn[j], conn2, count2, true );
            if( MB_SUCCESS != rval ) return rval;
 
            for( int k = 0; k < count2; ++k )
            {
                rval = moab()->get_coords( conn2 + k, 1, coords[0].array() );
                if( MB_SUCCESS != rval ) return rval;
                if( coords[0][plane.norm] <= plane.coord ) lo = true;
                if( coords[0][plane.norm] >= plane.coord ) ro = true;
            }
        }
 
        if( lo && ro )
            both_ins = both_tris.insert( both_ins, *i, *i );
        else if( lo )
            left_ins = left_tris.insert( left_ins, *i, *i );
        else if( ro )
            right_ins = right_tris.insert( right_ins, *i, *i );
    }
 
    // sets
    BoundBox tbox;
    for( i = set_begin; i != elems.end(); ++i )
    {
        tree_stats().constructLeafObjectTests++;
        rval = tbox.update( *moab(), *i, spherical, radius );
        if( MB_SUCCESS != rval ) return rval;
 
        bool lo = false, ro = false;
        if( tbox.bMin[plane.norm] <= plane.coord ) lo = true;
        if( tbox.bMax[plane.norm] >= plane.coord ) ro = true;
 
        if( lo && ro )
            both_ins = both_tris.insert( both_ins, *i, *i );
        else if( lo )
            left_ins = left_tris.insert( left_ins, *i, *i );
        else  // if (ro)
            right_ins = right_tris.insert( right_ins, *i, *i );
    }
 
    CartVect box_dim  = box_max - box_min;
    double area_left  = left_dim[0] * left_dim[1] + left_dim[1] * left_dim[2] + left_dim[2] * left_dim[0];
    double area_right = right_dim[0] * right_dim[1] + right_dim[1] * right_dim[2] + right_dim[2] * right_dim[0];
    double area_both  = box_dim[0] * box_dim[1] + box_dim[1] * box_dim[2] + box_dim[2] * box_dim[0];
    metric_value = ( area_left * left_tris.size() + area_right * right_tris.size() ) / area_both + both_tris.size();
    return MB_SUCCESS;
}

References moab::Range::begin(), moab::BoundBox::bMax, moab::BoundBox::bMin, moab::GeomUtil::box_elem_overlap(), box_max(), box_min(), moab::TreeStats::boxElemTests, moab::Range::clear(), moab::TreeStats::constructLeafObjectTests, moab::AdaptiveKDTree::Plane::coord, moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Range::insert(), moab::Range::lower_bound(), MB_SUCCESS, MBEDGE, MBENTITYSET, MBPOLYHEDRON, moab::Tree::moab(), moab::AdaptiveKDTree::Plane::norm, radius, moab::Range::size(), spherical, moab::Tree::tree_stats(), moab::TYPE_FROM_HANDLE(), and moab::BoundBox::update().

Referenced by best_subdivision_plane(), best_subdivision_snap_plane(), best_vertex_median_plane(), and best_vertex_sample_plane().

make_tag()

ErrorCode moab::AdaptiveKDTree::make_tag ( Interface *  iface, std::string  name, TagType  storage, DataType  type, int  count, void *  default_val, Tag &  tag_handle, std::vector< Tag > &  created_tags  )

private

Definition at line 209 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval =
        iface->tag_get_handle( name.c_str(), count, type, tag_handle, MB_TAG_CREAT | storage, default_val );
 
    if( MB_SUCCESS == rval )
    {
        if( std::find( created_tags.begin(), created_tags.end(), tag_handle ) == created_tags.end() )
            created_tags.push_back( tag_handle );
    }
    else
    {
        while( !created_tags.empty() )
        {
            iface->tag_delete( created_tags.back() );
            created_tags.pop_back();
        }
 
        planeTag = axisTag = (Tag)-1;
    }
 
    return rval;
}

References axisTag, ErrorCode, iface, MB_SUCCESS, MB_TAG_CREAT, and planeTag.

Referenced by init().

merge_leaf()

ErrorCode moab::AdaptiveKDTree::merge_leaf

(

AdaptiveKDTreeIter &

iter

)

Merge the leaf pointed to by the current iterator with it's sibling. If the sibling is not a leaf, multiple merges may be done.

Definition at line 419 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval;
    if( iter.depth() == 1 )  // at root
        return MB_FAILURE;
 
    // Move iter to parent
 
    AdaptiveKDTreeIter::StackObj node = iter.mStack.back();
    iter.mStack.pop_back();
 
    iter.childVect.clear();
    rval = moab()->get_child_meshsets( iter.mStack.back().entity, iter.childVect );
    if( MB_SUCCESS != rval ) return rval;
    Plane plane;
    rval = get_split_plane( iter.mStack.back().entity, plane );
    if( MB_SUCCESS != rval ) return rval;
 
    int child_idx = iter.childVect[0] == node.entity ? 0 : 1;
    assert( iter.childVect[child_idx] == node.entity );
    iter.mBox[1 - child_idx][plane.norm] = node.coord;
 
    // Get all entities from children and put them in parent
    EntityHandle parent = iter.handle();
    moab()->remove_child_meshset( parent, iter.childVect[0] );
    moab()->remove_child_meshset( parent, iter.childVect[1] );
    std::vector< EntityHandle > stack( iter.childVect );
 
    Range range;
    while( !stack.empty() )
    {
        EntityHandle h = stack.back();
        stack.pop_back();
        range.clear();
        rval = moab()->get_entities_by_handle( h, range );
        if( MB_SUCCESS != rval ) return rval;
        rval = moab()->add_entities( parent, range );
        if( MB_SUCCESS != rval ) return rval;
 
        iter.childVect.clear();
        MB_CHK_ERR( moab()->get_child_meshsets( h, iter.childVect ) );
        if( !iter.childVect.empty() )
        {
            moab()->remove_child_meshset( h, iter.childVect[0] );
            moab()->remove_child_meshset( h, iter.childVect[1] );
            stack.push_back( iter.childVect[0] );
            stack.push_back( iter.childVect[1] );
        }
 
        rval = moab()->delete_entities( &h, 1 );
        if( MB_SUCCESS != rval ) return rval;
    }
 
    return MB_SUCCESS;
}

References moab::Interface::add_entities(), moab::Range::back(), moab::AdaptiveKDTreeIter::childVect, moab::Range::clear(), moab::AdaptiveKDTreeIter::StackObj::coord, moab::Interface::delete_entities(), moab::AdaptiveKDTreeIter::depth(), moab::AdaptiveKDTreeIter::StackObj::entity, ErrorCode, moab::Interface::get_child_meshsets(), moab::Interface::get_entities_by_handle(), get_split_plane(), moab::AdaptiveKDTreeIter::handle(), MB_CHK_ERR, MB_SUCCESS, moab::AdaptiveKDTreeIter::mBox, moab::Tree::moab(), moab::AdaptiveKDTreeIter::mStack, moab::AdaptiveKDTree::Plane::norm, and moab::Interface::remove_child_meshset().

parse_options()

ErrorCode moab::AdaptiveKDTree::parse_options ( FileOptions &  options )

virtual

Parse options for tree creation.

Parameters

options

Options passed in by application

Returns

Failure is returned if any options were passed in and not interpreted; could mean inappropriate options for a particular tree type

Implements moab::Tree.

Definition at line 176 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval = parse_common_options( opts );
    if( MB_SUCCESS != rval ) return rval;
 
    //  SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
    int tmp_int;
    rval = opts.get_int_option( "SPLITS_PER_DIR", tmp_int );
    if( MB_SUCCESS == rval ) splitsPerDir = tmp_int;
 
    //  PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below)
    //           for possible values; default = 1 (SUBDIVISION_SNAP)
    rval = opts.get_int_option( "PLANE_SET", tmp_int );
    if( MB_SUCCESS == rval && ( tmp_int < SUBDIVISION || tmp_int > VERTEX_SAMPLE ) )
        return MB_FAILURE;
    else if( MB_ENTITY_NOT_FOUND == rval )
        planeSet = SUBDIVISION;
    else
        planeSet = (CandidatePlaneSet)( tmp_int );
 
    rval = opts.get_toggle_option( "SPHERICAL", false, spherical );
    if( MB_SUCCESS != rval ) spherical = false;
 
    double tmp = 1.0;
    rval       = opts.get_real_option( "RADIUS", tmp );
    if( MB_SUCCESS != rval )
        radius = 1.0;
    else
        radius = tmp;
 
    return MB_SUCCESS;
}

References ErrorCode, moab::FileOptions::get_int_option(), moab::FileOptions::get_real_option(), moab::FileOptions::get_toggle_option(), MB_ENTITY_NOT_FOUND, MB_SUCCESS, moab::Tree::parse_common_options(), planeSet, radius, spherical, splitsPerDir, SUBDIVISION, and VERTEX_SAMPLE.

Referenced by AdaptiveKDTree(), build_tree(), moab::Intx2Mesh::intersect_meshes(), and moab::Intx2Mesh::intersect_meshes_kdtree().

point_search() [1/2]

ErrorCode moab::AdaptiveKDTree::point_search	(	const double *	point,
		AdaptiveKDTreeIter &	leaf_it,
		const double	iter_tol = 1.0e-10,
		const double	inside_tol = 1.0e-6,
		bool *	multiple_leaves = NULL,
		EntityHandle *	start_node = NULL
	)

Get leaf containing input position.

Does not take into account global bounding box of tree.

• Therefore there is always one leaf containing the point.
• If caller wants to account for global bounding box, then caller can test against that box and not call this method at all if the point is outside the box, as there is no leaf containing the point in that case.

  Parameters

  point	Point to be located in tree
  leaf_it	Iterator to leaf containing point
  iter_tol	Tolerance for convergence of point search
  inside_tol	Tolerance for inside element calculation
  multiple_leaves	Some tree types can have multiple leaves containing a point; if non-NULL, this parameter is returned true if multiple leaves contain the input point
  start_node	Start from this tree node (non-NULL) instead of tree root (NULL)

  Returns

  Non-success returned only in case of failure; not-found indicated by leaf_out=0

Definition at line 1322 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval;
    treeStats.numTraversals++;
 
    // kdtrees never have multiple leaves containing a pt
    if( multiple_leaves ) *multiple_leaves = false;
 
    leaf_it.mBox[0] = boundBox.bMin;
    leaf_it.mBox[1] = boundBox.bMax;
 
    // test that point is inside tree
    if( !boundBox.contains_point( point, iter_tol ) )
    {
        treeStats.nodesVisited++;
        return MB_ENTITY_NOT_FOUND;
    }
 
    // initialize iterator at tree root
    leaf_it.treeTool = this;
    leaf_it.mStack.clear();
    leaf_it.mStack.push_back( AdaptiveKDTreeIter::StackObj( ( start_node ? *start_node : myRoot ), 0 ) );
 
    // loop until we reach a leaf
    AdaptiveKDTree::Plane plane;
    for( ;; )
    {
        treeStats.nodesVisited++;
 
        // get children
        leaf_it.childVect.clear();
        rval = moab()->get_child_meshsets( leaf_it.handle(), leaf_it.childVect );
        if( MB_SUCCESS != rval ) return rval;
 
        // if no children, then at leaf (done)
        if( leaf_it.childVect.empty() )
        {
            treeStats.leavesVisited++;
            break;
        }
 
        // get split plane
        rval = get_split_plane( leaf_it.handle(), plane );
        if( MB_SUCCESS != rval ) return rval;
 
        // step iterator to appropriate child
        // idx: 0->left, 1->right
        const int idx = ( point[plane.norm] > plane.coord );
        leaf_it.mStack.push_back(
            AdaptiveKDTreeIter::StackObj( leaf_it.childVect[idx], leaf_it.mBox[1 - idx][plane.norm] ) );
        leaf_it.mBox[1 - idx][plane.norm] = plane.coord;
    }
 
    return MB_SUCCESS;
}

References moab::BoundBox::bMax, moab::BoundBox::bMin, moab::Tree::boundBox, moab::AdaptiveKDTreeIter::childVect, moab::BoundBox::contains_point(), moab::AdaptiveKDTree::Plane::coord, ErrorCode, moab::Interface::get_child_meshsets(), get_split_plane(), moab::AdaptiveKDTreeIter::handle(), moab::TreeStats::leavesVisited, MB_ENTITY_NOT_FOUND, MB_SUCCESS, moab::AdaptiveKDTreeIter::mBox, moab::Tree::moab(), moab::AdaptiveKDTreeIter::mStack, moab::Tree::myRoot, moab::TreeStats::nodesVisited, moab::AdaptiveKDTree::Plane::norm, moab::TreeStats::numTraversals, moab::Tree::treeStats, and moab::AdaptiveKDTreeIter::treeTool.

point_search() [2/2]

ErrorCode moab::AdaptiveKDTree::point_search ( const double *  point, EntityHandle &  leaf_out, const double  iter_tol = 1.0e-10, const double  inside_tol = 1.0e-6, bool *  multiple_leaves = NULL, EntityHandle *  start_node = NULL, CartVect *  params = NULL  )

virtual

Get leaf containing input position.

Does not take into account global bounding box of tree.

• Therefore there is always one leaf containing the point.
• If caller wants to account for global bounding box, then caller can test against that box and not call this method at all if the point is outside the box, as there is no leaf containing the point in that case.

  Parameters

  point	Point to be located in tree
  leaf_out	Leaf containing point
  iter_tol	Tolerance for convergence of point search
  inside_tol	Tolerance for inside element calculation
  multiple_leaves	Some tree types can have multiple leaves containing a point; if non-NULL, this parameter is returned true if multiple leaves contain the input point
  start_node	Start from this tree node (non-NULL) instead of tree root (NULL)

  Returns

  Non-success returned only in case of failure; not-found indicated by leaf_out=0

Implements moab::Tree.

Definition at line 1267 of file AdaptiveKDTree.cpp.

{
    std::vector< EntityHandle > children;
    Plane plane;
 
    treeStats.numTraversals++;
    leaf_out = 0;
    BoundBox box;
    // kdtrees never have multiple leaves containing a pt
    if( multiple_leaves ) *multiple_leaves = false;
 
    EntityHandle node = ( start_node ? *start_node : myRoot );
 
    treeStats.nodesVisited++;
    ErrorCode rval = get_bounding_box( box, &node );
    if( MB_SUCCESS != rval ) return rval;
    if( !box.contains_point( point, iter_tol ) ) return MB_SUCCESS;
 
    rval = moab()->get_child_meshsets( node, children );
    if( MB_SUCCESS != rval ) return rval;
 
    while( !children.empty() )
    {
        treeStats.nodesVisited++;
 
        rval = get_split_plane( node, plane );
        if( MB_SUCCESS != rval ) return rval;
 
        const double d = point[plane.norm] - plane.coord;
        node           = children[( d > 0.0 )];
 
        children.clear();
        rval = moab()->get_child_meshsets( node, children );
        if( MB_SUCCESS != rval ) return rval;
    }
 
    treeStats.leavesVisited++;
    if( myEval && params )
    {
        rval = myEval->find_containing_entity( node, point, iter_tol, inside_tol, leaf_out, params->array(),
                                               &treeStats.traversalLeafObjectTests );
        if( MB_SUCCESS != rval ) return rval;
    }
    else
        leaf_out = node;
 
    return MB_SUCCESS;
}

References moab::BoundBox::contains_point(), moab::AdaptiveKDTree::Plane::coord, ErrorCode, moab::ElemEvaluator::find_containing_entity(), moab::Tree::get_bounding_box(), moab::Interface::get_child_meshsets(), get_split_plane(), moab::TreeStats::leavesVisited, MB_SUCCESS, moab::Tree::moab(), moab::Tree::myEval, moab::Tree::myRoot, moab::TreeStats::nodesVisited, moab::AdaptiveKDTree::Plane::norm, moab::TreeStats::numTraversals, moab::params, moab::TreeStats::traversalLeafObjectTests, and moab::Tree::treeStats.

Referenced by moab::TempestRemapper::ConstructCoveringSet(), main(), and moab::Coupler::nat_param().

print()

ErrorCode moab::AdaptiveKDTree::print ( )

virtual

print various things about this tree

Implements moab::Tree.

Definition at line 2130 of file AdaptiveKDTree.cpp.

{
    Range range;
 
    Range tree_sets, elem2d, elem3d, verts, all;
    moab()->get_child_meshsets( myRoot, tree_sets, 0 );
 
    {
        std::vector< EntityHandle > elem2d_vec, elem3d_vec, verts_vec;
        for( Range::iterator rit = tree_sets.begin(); rit != tree_sets.end(); ++rit )
        {
            moab()->get_entities_by_dimension( *rit, 2, elem2d_vec );
            moab()->get_entities_by_dimension( *rit, 3, elem3d_vec );
            moab()->get_entities_by_type( *rit, MBVERTEX, verts_vec );
        }
        std::sort( elem2d_vec.begin(), elem2d_vec.end() );
        std::copy( elem2d_vec.rbegin(), elem2d_vec.rend(), range_inserter( elem2d ) );
        std::sort( elem3d_vec.begin(), elem3d_vec.end() );
        std::copy( elem3d_vec.rbegin(), elem3d_vec.rend(), range_inserter( elem3d ) );
        std::sort( verts_vec.begin(), verts_vec.end() );
        std::copy( verts_vec.rbegin(), verts_vec.rend(), range_inserter( verts ) );
    }
 
    all.merge( verts );
    all.merge( elem2d );
    all.merge( elem3d );
    tree_sets.insert( myRoot );
    unsigned long long set_used, set_amortized, set_store_used, set_store_amortized, set_tag_used, set_tag_amortized,
        elem_used, elem_amortized;
    moab()->estimated_memory_use( tree_sets, &set_used, &set_amortized, &set_store_used, &set_store_amortized, 0, 0, 0,
                                  0, &set_tag_used, &set_tag_amortized );
    moab()->estimated_memory_use( all, &elem_used, &elem_amortized );
 
    int num_2d = 0, num_3d = 0;
    ;
    moab()->get_number_entities_by_dimension( 0, 2, num_2d );
    moab()->get_number_entities_by_dimension( 0, 3, num_3d );
 
    BoundBox box;
    ErrorCode rval = get_bounding_box( box, &myRoot );
    if( MB_SUCCESS != rval || box == BoundBox() ) throw rval;
    double diff[3]        = { box.bMax[0] - box.bMin[0], box.bMax[1] - box.bMin[1], box.bMax[2] - box.bMin[2] };
    double tree_vol       = diff[0] * diff[1] * diff[2];
    double tree_surf_area = 2 * ( diff[0] * diff[1] + diff[1] * diff[2] + diff[2] * diff[0] );
 
    SimpleStat< unsigned > depth, size;
    SimpleStat< double > vol, surf;
 
    AdaptiveKDTreeIter iter;
    get_tree_iterator( myRoot, iter );
    do
    {
        depth.add( iter.depth() );
 
        int num_leaf_elem;
        moab()->get_number_entities_by_handle( iter.handle(), num_leaf_elem );
        size.add( num_leaf_elem );
 
        const double* n = iter.box_min();
        const double* x = iter.box_max();
        double dims[3]  = { x[0] - n[0], x[1] - n[1], x[2] - n[2] };
 
        double leaf_vol = dims[0] * dims[1] * dims[2];
        vol.add( leaf_vol );
 
        double area = 2.0 * ( dims[0] * dims[1] + dims[1] * dims[2] + dims[2] * dims[0] );
        surf.add( area );
 
    } while( MB_SUCCESS == iter.step() );
 
    printf( "------------------------------------------------------------------\n" );
    printf( "tree volume:      %f\n", tree_vol );
    printf( "total elements:   %d\n", num_2d + num_3d );
    printf( "number of leaves: %lu\n", (unsigned long)depth.count );
    printf( "number of nodes:  %lu\n", (unsigned long)tree_sets.size() );
    printf( "volume ratio:     %0.2f%%\n", 100 * ( vol.sum / tree_vol ) );
    printf( "surface ratio:    %0.2f%%\n", 100 * ( surf.sum / tree_surf_area ) );
    printf( "\nmemory:           used  amortized\n" );
    printf( "            ---------- ----------\n" );
    printf( "elements    %10s %10s\n", mem_to_string( elem_used ).c_str(), mem_to_string( elem_amortized ).c_str() );
    printf( "sets (total)%10s %10s\n", mem_to_string( set_used ).c_str(), mem_to_string( set_amortized ).c_str() );
    printf( "sets        %10s %10s\n", mem_to_string( set_store_used ).c_str(),
            mem_to_string( set_store_amortized ).c_str() );
    printf( "set tags    %10s %10s\n", mem_to_string( set_tag_used ).c_str(),
            mem_to_string( set_tag_amortized ).c_str() );
    printf( "\nleaf stats:        min        avg        rms        max    std.dev\n" );
    printf( "            ---------- ---------- ---------- ---------- ----------\n" );
    printf( "depth       %10u %10.1f %10.1f %10u %10.2f\n", depth.min, depth.avg(), depth.rms(), depth.max,
            depth.dev() );
    printf( "triangles   %10u %10.1f %10.1f %10u %10.2f\n", size.min, size.avg(), size.rms(), size.max, size.dev() );
    printf( "volume      %10.2g %10.2g %10.2g %10.2g %10.2g\n", vol.min, vol.avg(), vol.rms(), vol.max, vol.dev() );
    printf( "surf. area  %10.2g %10.2g %10.2g %10.2g %10.2g\n", surf.min, surf.avg(), surf.rms(), surf.max,
            surf.dev() );
    printf( "------------------------------------------------------------------\n" );
 
    return MB_SUCCESS;
}

References moab::SimpleStat< T >::add(), moab::SimpleStat< T >::avg(), moab::Range::begin(), moab::BoundBox::bMax, moab::BoundBox::bMin, moab::AdaptiveKDTreeIter::box_max(), moab::AdaptiveKDTreeIter::box_min(), moab::SimpleStat< T >::count, moab::AdaptiveKDTreeIter::depth(), moab::SimpleStat< T >::dev(), moab::Range::end(), ErrorCode, moab::Interface::estimated_memory_use(), moab::Tree::get_bounding_box(), moab::Interface::get_child_meshsets(), moab::Interface::get_entities_by_dimension(), moab::Interface::get_entities_by_type(), moab::Interface::get_number_entities_by_dimension(), moab::Interface::get_number_entities_by_handle(), get_tree_iterator(), moab::AdaptiveKDTreeIter::handle(), moab::Range::insert(), moab::SimpleStat< T >::max, MB_SUCCESS, MBVERTEX, moab::mem_to_string(), moab::Range::merge(), moab::SimpleStat< T >::min, moab::Tree::moab(), moab::Tree::myRoot, moab::SimpleStat< T >::rms(), moab::Range::size(), moab::AdaptiveKDTreeIter::step(), and moab::SimpleStat< T >::sum.

ray_intersect_triangles()

ErrorCode moab::AdaptiveKDTree::ray_intersect_triangles	(	EntityHandle	tree_root,
		const double	tolerance,
		const double	ray_unit_dir[3],
		const double	ray_base_pt[3],
		std::vector< EntityHandle > &	triangles_out,
		std::vector< double > &	distance_out,
		int	result_count_limit = 0,
		double	distance_limit = -1.0
	)

Definition at line 1852 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval;
    double ray_beg = 0.0;
    if( ray_end < 0.0 ) ray_end = HUGE_VAL;
 
    // if root has bounding box, trim ray to that box
    CartVect tvec( tol );
    BoundBox box;
    const CartVect ray_pt( ray_pt_in ), ray_dir( ray_dir_in );
    rval = get_bounding_box( box );
    if( MB_SUCCESS == rval )
    {
        if( !GeomUtil::segment_box_intersect( box.bMin - tvec, box.bMax + tvec, ray_pt, ray_dir, ray_beg, ray_end ) )
            return MB_SUCCESS;  // ray misses entire tree.
    }
 
    Range tris;
    Range::iterator iter;
    CartVect tri_coords[3];
    const EntityHandle* tri_conn;
    int conn_len;
    double tri_t;
 
    Plane plane;
    std::vector< EntityHandle > children;
    std::vector< NodeSeg > list;
    NodeSeg seg( root, ray_beg, ray_end );
    list.push_back( seg );
 
    while( !list.empty() )
    {
        seg = list.back();
        list.pop_back();
 
        // If we are limited to a certain number of intersections
        // (max_ints != 0), then ray_end will contain the distance
        // to the furthest intersection we have so far.  If the
        // tree node is further than that, skip it.
        if( seg.beg > ray_end ) continue;
 
        // Check if at a leaf
        children.clear();
        rval = moab()->get_child_meshsets( seg.handle, children );
        if( MB_SUCCESS != rval ) return rval;
        if( children.empty() )
        {  // leaf
 
            tris.clear();
            rval = moab()->get_entities_by_type( seg.handle, MBTRI, tris );
            if( MB_SUCCESS != rval ) return rval;
 
            for( iter = tris.begin(); iter != tris.end(); ++iter )
            {
                rval = moab()->get_connectivity( *iter, tri_conn, conn_len );
                if( MB_SUCCESS != rval ) return rval;
                rval = moab()->get_coords( tri_conn, 3, tri_coords[0].array() );
                if( MB_SUCCESS != rval ) return rval;
 
                if( GeomUtil::ray_tri_intersect( tri_coords, ray_pt, ray_dir, tri_t, &ray_end ) )
                {
                    if( !max_ints )
                    {
                        if( std::find( tris_out.begin(), tris_out.end(), *iter ) == tris_out.end() )
                        {
                            tris_out.push_back( *iter );
                            dists_out.push_back( tri_t );
                        }
                    }
                    else if( tri_t < ray_end )
                    {
                        if( std::find( tris_out.begin(), tris_out.end(), *iter ) == tris_out.end() )
                        {
                            if( tris_out.size() < (unsigned)max_ints )
                            {
                                tris_out.resize( tris_out.size() + 1 );
                                dists_out.resize( dists_out.size() + 1 );
                            }
                            int w = tris_out.size() - 1;
                            for( ; w > 0 && tri_t < dists_out[w - 1]; --w )
                            {
                                tris_out[w]  = tris_out[w - 1];
                                dists_out[w] = dists_out[w - 1];
                            }
                            tris_out[w]  = *iter;
                            dists_out[w] = tri_t;
                            if( tris_out.size() >= (unsigned)max_ints )
                                // when we have already reached the max intx points, we cans safely
                                // reset ray_end, because we will accept new points only "closer"
                                // than the last one
                                ray_end = dists_out.back();
                        }
                    }
                }
            }
 
            continue;
        }
 
        rval = get_split_plane( seg.handle, plane );
        if( MB_SUCCESS != rval ) return rval;
 
        // Consider two planes that are the split plane +/- the tolerance.
        // Calculate the segment parameter at which the line segment intersects
        // the true plane, and also the difference between that value and the
        // intersection with either of the +/- tol planes.
        const double inv_dir = 1.0 / ray_dir[plane.norm];                       // only do division once
        const double t       = ( plane.coord - ray_pt[plane.norm] ) * inv_dir;  // intersection with plane
        const double diff    = tol * inv_dir;                                   // t adjustment for +tol plane
        // const double t0 = t - diff; // intersection with -tol plane
        // const double t1 = t + diff; // intersection with +tol plane
 
        // The index of the child tree node (0 or 1) that is on the
        // side of the plane to which the ray direction points.  That is,
        // if the ray direction is opposite the plane normal, the index
        // of the child corresponding to the side beneath the plane.  If
        // the ray direction is the same as the plane normal, the index
        // of the child corresponding to the side above the plane.
        const int fwd_child = ( ray_dir[plane.norm] > 0.0 );
 
        // Note: we maintain seg.beg <= seg.end at all times, so assume that here.
 
        // If segment is parallel to plane
        if( !Util::is_finite( t ) )
        {
            if( ray_pt[plane.norm] - tol <= plane.coord ) list.push_back( NodeSeg( children[0], seg.beg, seg.end ) );
            if( ray_pt[plane.norm] + tol >= plane.coord ) list.push_back( NodeSeg( children[1], seg.beg, seg.end ) );
        }
        // If segment is entirely to one side of plane such that the
        // intersection with the split plane is past the end of the segment
        else if( seg.end + diff < t )
        {
            // If segment direction is opposite that of plane normal, then
            // being past the end of the segment means that we are to the
            // right (or above) the plane and what the right child (index == 1).
            // Otherwise we want the left child (index == 0);
            list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.end ) );
        }
        // If the segment is entirely to one side of the plane such that
        // the intersection with the split plane is before the start of the
        // segment
        else if( seg.beg - diff > t )
        {
            // If segment direction is opposite that of plane normal, then
            // being before the start of the segment means that we are to the
            // left (or below) the plane and what the left child (index == 0).
            // Otherwise we want the right child (index == 1);
            list.push_back( NodeSeg( children[fwd_child], seg.beg, seg.end ) );
        }
        // Otherwise we must intersect the plane.
        // Note: be careful not to grow the segment if t is slightly
        // outside the current segment, as doing so would effectively
        // increase the tolerance as we descend the tree.
        else if( t <= seg.beg )
        {
            list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.beg ) );
            list.push_back( NodeSeg( children[fwd_child], seg.beg, seg.end ) );
        }
        else if( t >= seg.end )
        {
            list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.end ) );
            list.push_back( NodeSeg( children[fwd_child], seg.end, seg.end ) );
        }
        else
        {
            list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, t ) );
            list.push_back( NodeSeg( children[fwd_child], t, seg.end ) );
        }
    }
 
    return MB_SUCCESS;
}

References moab::NodeSeg::beg, moab::Range::begin(), moab::BoundBox::bMax, moab::BoundBox::bMin, moab::Range::clear(), moab::AdaptiveKDTree::Plane::coord, moab::NodeSeg::end, moab::Range::end(), ErrorCode, moab::Tree::get_bounding_box(), moab::Interface::get_child_meshsets(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), get_split_plane(), moab::NodeSeg::handle, moab::Util::is_finite(), MB_SUCCESS, MBTRI, moab::Tree::moab(), moab::AdaptiveKDTree::Plane::norm, moab::GeomUtil::ray_tri_intersect(), and moab::GeomUtil::segment_box_intersect().

reset_tree()

ErrorCode moab::AdaptiveKDTree::reset_tree ( )

inlinevirtual

Reset the tree, optionally checking we have the right root.

Implements moab::Tree.

Definition at line 528 of file AdaptiveKDTree.hpp.

{
    return delete_tree_sets();
}

References moab::Tree::delete_tree_sets().

Referenced by build_tree(), main(), and ~AdaptiveKDTree().

set_split_plane()

ErrorCode moab::AdaptiveKDTree::set_split_plane	(	EntityHandle	node,
		const Plane &	plane
	)

Set split plane for tree node.

Definition at line 303 of file AdaptiveKDTree.cpp.

{
#ifndef MB_AD_KD_TREE_USE_SINGLE_TAG
    ErrorCode r1, r2;
    r1 = moab()->tag_set_data( planeTag, &entity, 1, &plane.coord );
    r2 = moab()->tag_set_data( axisTag, &entity, 1, &plane.norm );
    return MB_SUCCESS == r1 ? r2 : r1;
#elif defined( MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG )
    double values[2] = { plane.coord, static_cast< double >( plane.norm ) };
    return moab()->tag_set_data( planeTag, &entity, 1, values );
#else
    return moab()->tag_set_data( planeTag, &entity, 1, &plane );
#endif
}

References axisTag, moab::AdaptiveKDTree::Plane::coord, ErrorCode, MB_SUCCESS, moab::Tree::moab(), moab::AdaptiveKDTree::Plane::norm, planeTag, and moab::Interface::tag_set_data().

Referenced by split_leaf().

sphere_intersect_triangles()

ErrorCode moab::AdaptiveKDTree::sphere_intersect_triangles	(	EntityHandle	tree_root,
		const double	center[3],
		double	radius,
		std::vector< EntityHandle > &	triangles
	)

Definition at line 1802 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval;
    std::vector< EntityHandle > leaves;
    const CartVect from( center );
    CartVect closest_pt;
    const EntityHandle* conn;
    CartVect coords[3];
    int conn_len;
 
    // get leaves of tree that intersect sphere
    assert( tree_root );
    rval = distance_search( center, rad, leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root );
    if( MB_SUCCESS != rval ) return rval;
 
    // search each leaf for triangles intersecting sphere
    for( unsigned i = 0; i < leaves.size(); ++i )
    {
        Range tris;
        rval = moab()->get_entities_by_type( leaves[i], MBTRI, tris );
        if( MB_SUCCESS != rval ) return rval;
 
        for( Range::iterator j = tris.begin(); j != tris.end(); ++j )
        {
            rval = moab()->get_connectivity( *j, conn, conn_len );
            if( MB_SUCCESS != rval ) return rval;
            rval = moab()->get_coords( conn, 3, coords[0].array() );
            if( MB_SUCCESS != rval ) return rval;
            GeomUtil::closest_location_on_tri( from, coords, closest_pt );
            closest_pt -= from;
            if( ( closest_pt % closest_pt ) <= ( rad * rad ) ) triangles.push_back( *j );
        }
    }
 
    // remove duplicates from triangle list
    std::sort( triangles.begin(), triangles.end() );
    triangles.erase( std::unique( triangles.begin(), triangles.end() ), triangles.end() );
    return MB_SUCCESS;
}

References moab::Range::begin(), center(), moab::GeomUtil::closest_location_on_tri(), distance_search(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), MB_SUCCESS, MBTRI, and moab::Tree::moab().

split_leaf() [1/4]

ErrorCode moab::AdaptiveKDTree::split_leaf	(	AdaptiveKDTreeIter &	leaf,
		Plane	plane
	)

Split leaf of tree Updates iterator location to point to first new leaf node.

Definition at line 371 of file AdaptiveKDTree.cpp.

{
    EntityHandle left, right;
    return split_leaf( leaf, plane, left, right );
}

Referenced by build_tree(), and split_leaf().

split_leaf() [2/4]

ErrorCode moab::AdaptiveKDTree::split_leaf	(	AdaptiveKDTreeIter &	leaf,
		Plane	plane,
		const Range &	left_entities,
		const Range &	right_entities
	)

Split leaf of tree Updates iterator location to point to first new leaf node.

Definition at line 377 of file AdaptiveKDTree.cpp.

{
    EntityHandle left, right, parent = leaf.handle();
    ErrorCode rval = split_leaf( leaf, plane, left, right );
    if( MB_SUCCESS != rval ) return rval;
 
    if( MB_SUCCESS == moab()->add_entities( left, left_entities ) &&
        MB_SUCCESS == moab()->add_entities( right, right_entities ) &&
        MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
        return MB_SUCCESS;
 
    moab()->remove_child_meshset( parent, left );
    moab()->remove_child_meshset( parent, right );
    EntityHandle children[] = { left, right };
    moab()->delete_entities( children, 2 );
    return MB_FAILURE;
}

References moab::Interface::delete_entities(), ErrorCode, moab::AdaptiveKDTreeIter::handle(), MB_SUCCESS, moab::Tree::moab(), moab::Interface::remove_child_meshset(), and split_leaf().

split_leaf() [3/4]

ErrorCode moab::AdaptiveKDTree::split_leaf	(	AdaptiveKDTreeIter &	leaf,
		Plane	plane,
		const std::vector< EntityHandle > &	left_entities,
		const std::vector< EntityHandle > &	right_entities
	)

Split leaf of tree Updates iterator location to point to first new leaf node.

Definition at line 398 of file AdaptiveKDTree.cpp.

{
    EntityHandle left, right, parent = leaf.handle();
    ErrorCode rval = split_leaf( leaf, plane, left, right );
    if( MB_SUCCESS != rval ) return rval;
 
    if( MB_SUCCESS == moab()->add_entities( left, &left_entities[0], left_entities.size() ) &&
        MB_SUCCESS == moab()->add_entities( right, &right_entities[0], right_entities.size() ) &&
        MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
        return MB_SUCCESS;
 
    moab()->remove_child_meshset( parent, left );
    moab()->remove_child_meshset( parent, right );
    EntityHandle children[] = { left, right };
    moab()->delete_entities( children, 2 );
    return MB_FAILURE;
}

References moab::Interface::delete_entities(), ErrorCode, moab::AdaptiveKDTreeIter::handle(), MB_SUCCESS, moab::Tree::moab(), moab::Interface::remove_child_meshset(), and split_leaf().

split_leaf() [4/4]

ErrorCode moab::AdaptiveKDTree::split_leaf	(	AdaptiveKDTreeIter &	leaf,
		Plane	plane,
		EntityHandle &	left_child,
		EntityHandle &	right_child
	)

Split leaf of tree Updates iterator location to point to first new leaf node.

Definition at line 344 of file AdaptiveKDTree.cpp.

{
    ErrorCode rval;
 
    rval = moab()->create_meshset( meshsetFlags, left );
    if( MB_SUCCESS != rval ) return rval;
 
    rval = moab()->create_meshset( meshsetFlags, right );
    if( MB_SUCCESS != rval )
    {
        moab()->delete_entities( &left, 1 );
        return rval;
    }
 
    if( MB_SUCCESS != set_split_plane( leaf.handle(), plane ) ||
        MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), left ) ||
        MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), right ) ||
        MB_SUCCESS != leaf.step_to_first_leaf( AdaptiveKDTreeIter::LEFT ) )
    {
        EntityHandle children[] = { left, right };
        moab()->delete_entities( children, 2 );
        return MB_FAILURE;
    }
 
    return MB_SUCCESS;
}

References moab::Interface::create_meshset(), moab::Interface::delete_entities(), ErrorCode, moab::AdaptiveKDTreeIter::handle(), moab::AdaptiveKDTreeIter::LEFT, MB_SUCCESS, moab::Tree::meshsetFlags, moab::Tree::moab(), set_split_plane(), and moab::AdaptiveKDTreeIter::step_to_first_leaf().

Friends And Related Function Documentation

AdaptiveKDTreeIter

friend class AdaptiveKDTreeIter

friend

Definition at line 264 of file AdaptiveKDTree.hpp.

Member Data Documentation

axisTag

Tag moab::AdaptiveKDTree::axisTag

private

Definition at line 331 of file AdaptiveKDTree.hpp.

Referenced by get_split_plane(), init(), make_tag(), and set_split_plane().

planeSet

CandidatePlaneSet moab::AdaptiveKDTree::planeSet

private

Definition at line 335 of file AdaptiveKDTree.hpp.

Referenced by build_tree(), and parse_options().

planeTag

Tag moab::AdaptiveKDTree::planeTag

private

Definition at line 331 of file AdaptiveKDTree.hpp.

Referenced by get_split_plane(), init(), make_tag(), and set_split_plane().

radius

double moab::AdaptiveKDTree::radius

private

Definition at line 338 of file AdaptiveKDTree.hpp.

Referenced by build_tree(), intersect_children_with_elems(), and parse_options().

spherical

bool moab::AdaptiveKDTree::spherical

private

Definition at line 337 of file AdaptiveKDTree.hpp.

Referenced by build_tree(), intersect_children_with_elems(), and parse_options().

splitsPerDir

unsigned moab::AdaptiveKDTree::splitsPerDir

private

Definition at line 333 of file AdaptiveKDTree.hpp.

Referenced by build_tree(), and parse_options().

treeName

const char * moab::AdaptiveKDTree::treeName = "AKDTree"

staticprivate

Definition at line 329 of file AdaptiveKDTree.hpp.

Referenced by AdaptiveKDTree(), and init().

---

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

• AdaptiveKDTree.hpp
• AdaptiveKDTree.cpp