point_location.cpp File Reference

#include "moab/Core.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/Range.hpp"
#include "moab/CartVect.hpp"
#include "moab/GeomUtil.hpp"
#include <iostream>
#include <ctime>
#include <cstdlib>
#include <cassert>
#include <sstream>

Include dependency graph for point_location.cpp:

Go to the source code of this file.

Macros
#define	CHK(ErrorCode)

Enumerations
enum	TreeType { UseKDTree , UseNoTree , UseDefaultTree = UseKDTree }

Functions
static void	fail (ErrorCode error_code, const char *file_name, int line_number)
const char *	is_default_tree (TreeType type)
static void	usage (char *argv0, bool help=false)
static void	generate_random_points (Interface &mesh, size_t num_points, std::vector< CartVect > &points, std::vector< EntityHandle > &point_elems)
static void	do_kdtree_test (Interface &mesh, int tree_depth, int elem_per_leaf, long num_test, const std::vector< CartVect > &points, std::vector< EntityHandle > &point_elems, clock_t &build_time, clock_t &test_time, size_t &depth)
static void	do_linear_test (Interface &mesh, int tree_depth, int elem_per_leaf, long num_test, const std::vector< CartVect > &points, std::vector< EntityHandle > &point_elems, clock_t &build_time, clock_t &test_time, size_t &depth)
int	main (int argc, char *argv[])
static CartVect	random_point_in_hex (Interface &mb, EntityHandle hex)

Variables
const char *	default_str = "(default)"
const char *	empty_str = ""
const long	DEFAULT_NUM_TEST = 100000
const long	HARD_MAX_UNIQUE_POINTS = 100000
const long	HARD_MIN_UNIQUE_POINTS = 1000
const long	FRACTION_UNIQUE_POINTS = 100
const int	HexSign [8][3]

Macro Definition Documentation

CHK

#define CHK

(

ErrorCode

)

Value:

 do \
 { \
 if( MB_SUCCESS != ( ErrorCode ) ) fail( ( ErrorCode ), __FILE__, __LINE__ ); \
 } while( false )

Definition at line 12 of file point_location.cpp.

Enumeration Type Documentation

TreeType

enum TreeType

Enumerator
UseKDTree
UseNoTree
UseDefaultTree

Definition at line 22 of file point_location.cpp.

{
 UseKDTree,
 UseNoTree,
 UseDefaultTree = UseKDTree
};

Function Documentation

do_kdtree_test()

void do_kdtree_test ( Interface &  mesh, int  tree_depth, int  elem_per_leaf, long  num_test, const std::vector< CartVect > &  points, std::vector< EntityHandle > &  point_elems, clock_t &  build_time, clock_t &  test_time, size_t &  depth  )

static

Definition at line 310 of file point_location.cpp.

{
 ErrorCode rval;
 clock_t init = clock();
 AdaptiveKDTree tool( &mb );
 EntityHandle root;
 std::ostringstream options;
 if( tree_depth > 0 ) options << "MAX_DEPTH=" << tree_depth << ";";
 if( elem_per_leaf > 0 ) options << "MAX_PER_LEAF=" << elem_per_leaf << ";";
 Range all_hexes;
 rval = mb.get_entities_by_type( 0, MBHEX, all_hexes );
 CHK( rval );
 FileOptions opt( options.str().c_str() );
 rval = tool.build_tree( all_hexes, &root, &opt );
 CHK( rval );
 all_hexes.clear();
 build_time = clock() - init;
 
 EntityHandle leaf;
 std::vector< EntityHandle > hexes;
 std::vector< EntityHandle >::iterator j;
 CartVect coords[8];
 const EntityHandle* conn;
 int len;
 for( long i = 0; i < num_test; ++i )
 {
 const size_t idx = (size_t)i % points.size();
 rval = tool.point_search( points[idx].array(), leaf, 1.0e-10, 1.0e-6, NULL, &root );
 CHK( rval );
 hexes.clear();
 rval = mb.get_entities_by_handle( leaf, hexes );
 CHK( rval );
 for( j = hexes.begin(); j != hexes.end(); ++j )
 {
 rval = mb.get_connectivity( *j, conn, len, true );
 CHK( rval );
 rval = mb.get_coords( conn, 8, reinterpret_cast< double* >( coords ) );
 CHK( rval );
 if( GeomUtil::point_in_trilinear_hex( coords, points[idx], 1e-12 ) )
 {
 point_elems[idx] = *j;
 break;
 }
 }
 if( j == hexes.end() ) point_elems[idx] = 0;
 }
 
 test_time = clock() - build_time;
 double tmp_box[3];
 unsigned max_d;
 tool.get_info( root, tmp_box, tmp_box, max_d );
 depth = max_d;
}

References build_time, moab::AdaptiveKDTree::build_tree(), CHK, moab::Range::clear(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), moab::Core::get_entities_by_handle(), moab::Core::get_entities_by_type(), moab::AdaptiveKDTree::get_info(), init(), leaf, mb, MBHEX, moab::GeomUtil::point_in_trilinear_hex(), and moab::AdaptiveKDTree::point_search().

Referenced by main().

do_linear_test()

void do_linear_test ( Interface &  mesh, int  tree_depth, int  elem_per_leaf, long  num_test, const std::vector< CartVect > &  points, std::vector< EntityHandle > &  point_elems, clock_t &  build_time, clock_t &  test_time, size_t &  depth  )

static

Definition at line 372 of file point_location.cpp.

{
 clock_t init = clock();
 Range hexes;
 ErrorCode rval = mb.get_entities_by_type( 0, MBHEX, hexes );
 CHK( rval );
 depth = 0;
 point_elems.resize( points.size() );
 build_time = clock() - init;
 
 CartVect coords[8];
 const EntityHandle* conn;
 int len;
 for( long i = 0; i < num_test; ++i )
 {
 const size_t idx = (size_t)i % points.size();
 for( Range::iterator h = hexes.begin(); h != hexes.end(); ++h )
 {
 rval = mb.get_connectivity( *h, conn, len, true );
 CHK( rval );
 rval = mb.get_coords( conn, 8, reinterpret_cast< double* >( coords ) );
 CHK( rval );
 if( GeomUtil::point_in_trilinear_hex( coords, points[idx], 1e-12 ) )
 {
 point_elems[idx] = *h;
 break;
 }
 }
 }
 
 test_time = clock() - build_time;
}

References moab::Range::begin(), build_time, CHK, moab::Range::end(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), moab::Core::get_entities_by_type(), init(), mb, MBHEX, and moab::GeomUtil::point_in_trilinear_hex().

Referenced by main().

fail()

void fail ( ErrorCode  error_code, const char *  file_name, int  line_number  )

static

Definition at line 236 of file point_location.cpp.

{
 std::cerr << "Internal error (error code " << error_code << ") at " << file << ":" << line << std::endl;
 abort();
}

generate_random_points()

void generate_random_points ( Interface &  mesh, size_t  num_points, std::vector< CartVect > &  points, std::vector< EntityHandle > &  point_elems  )

static

Definition at line 274 of file point_location.cpp.

{
 Range elems;
 ErrorCode rval;
 rval = mb.get_entities_by_dimension( 0, 3, elems );
 CHK( rval );
 if( !elems.all_of_type( MBHEX ) )
 {
 std::cerr << "Warning: ignoring non-hexahedral elements." << std::endl;
 std::pair< Range::iterator, Range::iterator > p = elems.equal_range( MBHEX );
 elems.erase( p.second, elems.end() );
 elems.erase( elems.begin(), p.first );
 }
 if( elems.empty() )
 {
 std::cerr << "Input file contains no hexahedral elements." << std::endl;
 exit( 1 );
 }
 
 points.resize( num_points );
 point_elems.resize( num_points );
 const size_t num_elem = elems.size();
 for( size_t i = 0; i < num_points; ++i )
 {
 size_t offset = 0;
 for( size_t x = num_elem; x > 0; x /= RAND_MAX )
 offset += rand();
 offset %= num_elem;
 point_elems[i] = elems[offset];
 points[i] = random_point_in_hex( mb, point_elems[i] );
 }
}

References moab::Range::all_of_type(), moab::Range::begin(), CHK, moab::Range::empty(), moab::Range::end(), moab::Range::equal_range(), moab::Range::erase(), ErrorCode, moab::Core::get_entities_by_dimension(), mb, MBHEX, random_point_in_hex(), and moab::Range::size().

Referenced by main().

is_default_tree()

const char* is_default_tree ( TreeType  type )

inline

Definition at line 31 of file point_location.cpp.

{
 return type == UseDefaultTree ? default_str : empty_str;
}

References default_str, empty_str, and UseDefaultTree.

Referenced by usage().

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 87 of file point_location.cpp.

{
 
 const char* input_file = 0;
 long tree_depth = -1;
 long elem_per_leaf = -1;
 long num_points = DEFAULT_NUM_TEST;
 TreeType type = UseDefaultTree;
 char* endptr;
 
 // PARSE ARGUMENTS
 
 long* valptr = 0;
 for( int i = 1; i < argc; ++i )
 {
 if( valptr )
 {
 *valptr = strtol( argv[i], &endptr, 0 );
 if( *valptr < 1 || *endptr )
 {
 std::cerr << "Invalid value for '" << argv[i - 1] << "' flag: " << argv[i] << std::endl;
 exit( 1 );
 }
 valptr = 0;
 }
 else if( argv[i][0] == '-' && argv[i][1] && !argv[i][2] )
 {
 switch( argv[i][1] )
 {
 case 'h':
 usage( argv[0], true );
 break;
 case 'k':
 type = UseKDTree;
 break;
 case 'v':
 type = UseNoTree;
 break;
 case 'd':
 valptr = &tree_depth;
 break;
 case 'e':
 valptr = &elem_per_leaf;
 break;
 case 'n':
 valptr = &num_points;
 break;
 default:
 std::cerr << "Invalid flag: " << argv[i] << std::endl;
 usage( argv[0] );
 break;
 }
 }
 else if( !input_file )
 {
 input_file = argv[i];
 }
 else
 {
 std::cerr << "Unexpected argument: " << argv[i] << std::endl;
 usage( argv[0] );
 }
 }
 if( valptr )
 {
 std::cerr << "Expected value following '" << argv[argc - 1] << "' flag" << std::endl;
 usage( argv[0] );
 }
 if( !input_file )
 {
 std::cerr << "No input file specified." << std::endl;
 usage( argv[0] );
 }
 
 // LOAD MESH
 
 Core moab;
 Interface& mb = moab;
 ErrorCode rval;
 std::string init_msg, msg;
 mb.get_last_error( init_msg );
 rval = mb.load_file( input_file );
 if( MB_SUCCESS != rval )
 {
 std::cerr << input_file << " : failed to read file." << std::endl;
 mb.get_last_error( msg );
 if( msg != init_msg ) std::cerr << "message : " << msg << std::endl;
 }
 
 // GENERATE TEST POINTS
 
 int num_unique = num_points / FRACTION_UNIQUE_POINTS;
 if( num_unique > HARD_MAX_UNIQUE_POINTS )
 num_unique = HARD_MAX_UNIQUE_POINTS;
 else if( num_unique < HARD_MIN_UNIQUE_POINTS )
 num_unique = num_points;
 std::vector< CartVect > points;
 std::vector< EntityHandle > elems;
 generate_random_points( mb, num_unique, points, elems );
 
 // GET MEMORY USE BEFORE BUILDING TREE
 
 unsigned long long init_total_storage;
 mb.estimated_memory_use( 0, 0, &init_total_storage );
 
 // RUN TIMING TEST
 clock_t build_time, test_time;
 size_t actual_depth = 0;
 std::vector< EntityHandle > results( points.size() );
 switch( type )
 {
 default:
 case UseKDTree:
 do_kdtree_test( mb, tree_depth, elem_per_leaf, num_points, points, results, build_time, test_time,
 actual_depth );
 break;
 case UseNoTree:
 do_linear_test( mb, tree_depth, elem_per_leaf, num_points, points, results, build_time, test_time,
 actual_depth );
 break;
 }
 
 unsigned long long fini_total_storage;
 mb.estimated_memory_use( 0, 0, &fini_total_storage );
 
 // VALIDATE RESULTS
 int fail = 0;
 if( results != elems )
 {
 std::cout << "WARNING: Test produced invalid results" << std::endl;
 fail = 1;
 }
 
 // SUMMARIZE RESULTS
 std::cout << "Number of test queries: " << num_points << std::endl;
 std::cout << "Tree build time: " << ( (double)build_time ) / CLOCKS_PER_SEC << " seconds" << std::endl;
 std::cout << "Total query time: " << ( (double)test_time ) / CLOCKS_PER_SEC << " seconds" << std::endl;
 std::cout << "Time per query: " << ( (double)test_time ) / CLOCKS_PER_SEC / num_points << " seconds" << std::endl;
 std::cout << "Tree depth: " << actual_depth << std::endl;
 if( actual_depth )
 std::cout << "Total query time/depth: " << ( (double)test_time ) / CLOCKS_PER_SEC / actual_depth << " seconds"
 << std::endl;
 std::cout << std::endl;
 std::cout << "Bytes before tree construction: " << init_total_storage << std::endl;
 std::cout << "Bytes after tree construction: " << fini_total_storage << std::endl;
 std::cout << "Difference: " << fini_total_storage - init_total_storage << " bytes" << std::endl;
 return fail;
}

References build_time, DEFAULT_NUM_TEST, do_kdtree_test(), do_linear_test(), ErrorCode, moab::Core::estimated_memory_use(), moab::fail(), FRACTION_UNIQUE_POINTS, generate_random_points(), moab::Core::get_last_error(), HARD_MAX_UNIQUE_POINTS, HARD_MIN_UNIQUE_POINTS, input_file, moab::Core::load_file(), mb, MB_SUCCESS, usage(), UseDefaultTree, UseKDTree, and UseNoTree.

random_point_in_hex()

static CartVect random_point_in_hex ( Interface &  mb, EntityHandle  hex  )

static

Definition at line 245 of file point_location.cpp.

{
 const double f = RAND_MAX / 2;
 CartVect xi( ( (double)rand() - f ) / f, ( (double)rand() - f ) / f, ( (double)rand() - f ) / f );
 CartVect coords[8];
 const EntityHandle* conn;
 int len;
 ErrorCode rval = mb.get_connectivity( hex, conn, len, true );
 if( len != 8 && MB_SUCCESS != rval )
 {
 std::cerr << "Invalid element" << std::endl;
 assert( false );
 abort();
 }
 rval = mb.get_coords( conn, 8, reinterpret_cast< double* >( coords ) );
 CHK( rval );
 
 CartVect point( 0, 0, 0 );
 for( unsigned i = 0; i < 8; ++i )
 {
 double coeff = 0.125;
 for( unsigned j = 0; j < 3; ++j )
 coeff *= 1 + HexSign[i][j] * xi[j];
 point += coeff * coords[i];
 }
 
 return point;
}

References CHK, ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), HexSign, mb, and MB_SUCCESS.

Referenced by generate_random_points().

usage()

static void usage ( char *  argv0, bool  help = false  )

static

Definition at line 41 of file point_location.cpp.

{
 const char* usage_str = "[-k|-v] [-n <N>] [-d <N>] [-e <N>] <input_mesh>";
 if( !help )
 {
 std::cerr << "Usage: " << argv0 << " " << usage_str << std::endl;
 std::cerr << " : " << argv0 << " -h" << std::endl;
 exit( 1 );
 }
 
 std::cout << "Usage: " << argv0 << " " << usage_str << std::endl
 << " -k : Use kD Tree " << is_default_tree( UseKDTree ) << std::endl
 << " -v : Use no tree" << is_default_tree( UseNoTree ) << std::endl
 << " -n : Specify number of test points (default = " << DEFAULT_NUM_TEST << ")" << std::endl
 << " -d : Specify maximum tree depth" << std::endl
 << " -e : Specify maximum elements per leaf" << std::endl
 << " <input_mesh> : Mesh to build and query." << std::endl
 << std::endl;
 exit( 0 );
}

References argv0, DEFAULT_NUM_TEST, help(), is_default_tree(), UseKDTree, and UseNoTree.

Referenced by main().

Variable Documentation

DEFAULT_NUM_TEST

const long DEFAULT_NUM_TEST = 100000

Definition at line 36 of file point_location.cpp.

Referenced by main(), and usage().

default_str

const char* default_str = "(default)"

Definition at line 29 of file point_location.cpp.

Referenced by is_default_tree().

empty_str

const char* empty_str = ""

Definition at line 30 of file point_location.cpp.

Referenced by is_default_tree().

FRACTION_UNIQUE_POINTS

const long FRACTION_UNIQUE_POINTS = 100

Definition at line 39 of file point_location.cpp.

Referenced by main().

HARD_MAX_UNIQUE_POINTS

const long HARD_MAX_UNIQUE_POINTS = 100000

Definition at line 37 of file point_location.cpp.

Referenced by main().

HARD_MIN_UNIQUE_POINTS

const long HARD_MIN_UNIQUE_POINTS = 1000

Definition at line 38 of file point_location.cpp.

Referenced by main().

HexSign

const int HexSign[8][3]

Initial value:

= { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
 { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } }

Definition at line 242 of file point_location.cpp.

Referenced by random_point_in_hex().