MOAB: Mesh Oriented datABase  (version 5.3.1)
point_location.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/AdaptiveKDTree.hpp"
00003 #include "moab/Range.hpp"
00004 #include "moab/CartVect.hpp"
00005 #include "moab/GeomUtil.hpp"
00006 #include <iostream>
00007 #include <ctime>
00008 #include <cstdlib>
00009 #include <cassert>
00010 #include <sstream>
00011 
00012 #define CHK( ErrorCode )                                                             \
00013     do                                                                               \
00014     {                                                                                \
00015         if( MB_SUCCESS != ( ErrorCode ) ) fail( ( ErrorCode ), __FILE__, __LINE__ ); \
00016     } while( false )
00017 
00018 using namespace moab;
00019 
00020 static void fail( ErrorCode error_code, const char* file_name, int line_number );
00021 
00022 enum TreeType
00023 {
00024     UseKDTree,
00025     UseNoTree,
00026     UseDefaultTree = UseKDTree
00027 };
00028 
00029 const char* default_str = "(default)";
00030 const char* empty_str   = "";
00031 inline const char* is_default_tree( TreeType type )
00032 {
00033     return type == UseDefaultTree ? default_str : empty_str;
00034 }
00035 
00036 const long DEFAULT_NUM_TEST       = 100000;
00037 const long HARD_MAX_UNIQUE_POINTS = 100000;
00038 const long HARD_MIN_UNIQUE_POINTS = 1000;
00039 const long FRACTION_UNIQUE_POINTS = 100;
00040 
00041 static void usage( char* argv0, bool help = false )
00042 {
00043     const char* usage_str = "[-k|-v] [-n <N>] [-d <N>] [-e <N>] <input_mesh>";
00044     if( !help )
00045     {
00046         std::cerr << "Usage: " << argv0 << " " << usage_str << std::endl;
00047         std::cerr << "     : " << argv0 << " -h" << std::endl;
00048         exit( 1 );
00049     }
00050 
00051     std::cout << "Usage: " << argv0 << " " << usage_str << std::endl
00052               << "  -k : Use kD Tree " << is_default_tree( UseKDTree ) << std::endl
00053               << "  -v : Use no tree" << is_default_tree( UseNoTree ) << std::endl
00054               << "  -n : Specify number of test points (default = " << DEFAULT_NUM_TEST << ")" << std::endl
00055               << "  -d : Specify maximum tree depth" << std::endl
00056               << "  -e : Specify maximum elements per leaf" << std::endl
00057               << " <input_mesh> : Mesh to build and query." << std::endl
00058               << std::endl;
00059     exit( 0 );
00060 }
00061 
00062 static void generate_random_points( Interface& mesh, size_t num_points, std::vector< CartVect >& points,
00063                                     std::vector< EntityHandle >& point_elems );
00064 
00065 static void do_kdtree_test( Interface& mesh, int tree_depth, int elem_per_leaf, long num_test,
00066                             const std::vector< CartVect >& points, std::vector< EntityHandle >& point_elems,
00067                             clock_t& build_time, clock_t& test_time, size_t& depth );
00068 
00069 static void do_linear_test( Interface& mesh, int tree_depth, int elem_per_leaf, long num_test,
00070                             const std::vector< CartVect >& points, std::vector< EntityHandle >& point_elems,
00071                             clock_t& build_time, clock_t& test_time, size_t& depth );
00072 
00073 int main( int argc, char* argv[] )
00074 {
00075 
00076     const char* input_file = 0;
00077     long tree_depth        = -1;
00078     long elem_per_leaf     = -1;
00079     long num_points        = DEFAULT_NUM_TEST;
00080     TreeType type          = UseDefaultTree;
00081     char* endptr;
00082 
00083     // PARSE ARGUMENTS
00084 
00085     long* valptr = 0;
00086     for( int i = 1; i < argc; ++i )
00087     {
00088         if( valptr )
00089         {
00090             *valptr = strtol( argv[i], &endptr, 0 );
00091             if( *valptr < 1 || *endptr )
00092             {
00093                 std::cerr << "Invalid value for '" << argv[i - 1] << "' flag: " << argv[i] << std::endl;
00094                 exit( 1 );
00095             }
00096             valptr = 0;
00097         }
00098         else if( argv[i][0] == '-' && argv[i][1] && !argv[i][2] )
00099         {
00100             switch( argv[i][1] )
00101             {
00102                 case 'h':
00103                     usage( argv[0], true );
00104                     break;
00105                 case 'k':
00106                     type = UseKDTree;
00107                     break;
00108                 case 'v':
00109                     type = UseNoTree;
00110                     break;
00111                 case 'd':
00112                     valptr = &tree_depth;
00113                     break;
00114                 case 'e':
00115                     valptr = &elem_per_leaf;
00116                     break;
00117                 case 'n':
00118                     valptr = &num_points;
00119                     break;
00120                 default:
00121                     std::cerr << "Invalid flag: " << argv[i] << std::endl;
00122                     usage( argv[0] );
00123                     break;
00124             }
00125         }
00126         else if( !input_file )
00127         {
00128             input_file = argv[i];
00129         }
00130         else
00131         {
00132             std::cerr << "Unexpected argument: " << argv[i] << std::endl;
00133             usage( argv[0] );
00134         }
00135     }
00136     if( valptr )
00137     {
00138         std::cerr << "Expected value following '" << argv[argc - 1] << "' flag" << std::endl;
00139         usage( argv[0] );
00140     }
00141     if( !input_file )
00142     {
00143         std::cerr << "No input file specified." << std::endl;
00144         usage( argv[0] );
00145     }
00146 
00147     // LOAD MESH
00148 
00149     Core moab;
00150     Interface& mb = moab;
00151     ErrorCode rval;
00152     std::string init_msg, msg;
00153     mb.get_last_error( init_msg );
00154     rval = mb.load_file( input_file );
00155     if( MB_SUCCESS != rval )
00156     {
00157         std::cerr << input_file << " : failed to read file." << std::endl;
00158         mb.get_last_error( msg );
00159         if( msg != init_msg ) std::cerr << "message : " << msg << std::endl;
00160     }
00161 
00162     // GENERATE TEST POINTS
00163 
00164     int num_unique = num_points / FRACTION_UNIQUE_POINTS;
00165     if( num_unique > HARD_MAX_UNIQUE_POINTS )
00166         num_unique = HARD_MAX_UNIQUE_POINTS;
00167     else if( num_unique < HARD_MIN_UNIQUE_POINTS )
00168         num_unique = num_points;
00169     std::vector< CartVect > points;
00170     std::vector< EntityHandle > elems;
00171     generate_random_points( mb, num_unique, points, elems );
00172 
00173     // GET MEMORY USE BEFORE BUILDING TREE
00174 
00175     unsigned long long init_total_storage;
00176     mb.estimated_memory_use( 0, 0, &init_total_storage );
00177 
00178     // RUN TIMING TEST
00179     clock_t build_time, test_time;
00180     size_t actual_depth = 0;
00181     std::vector< EntityHandle > results( points.size() );
00182     switch( type )
00183     {
00184         default:
00185         case UseKDTree:
00186             do_kdtree_test( mb, tree_depth, elem_per_leaf, num_points, points, results, build_time, test_time,
00187                             actual_depth );
00188             break;
00189         case UseNoTree:
00190             do_linear_test( mb, tree_depth, elem_per_leaf, num_points, points, results, build_time, test_time,
00191                             actual_depth );
00192             break;
00193     }
00194 
00195     unsigned long long fini_total_storage;
00196     mb.estimated_memory_use( 0, 0, &fini_total_storage );
00197 
00198     // VALIDATE RESULTS
00199     int fail = 0;
00200     if( results != elems )
00201     {
00202         std::cout << "WARNING:  Test produced invalid results" << std::endl;
00203         fail = 1;
00204     }
00205 
00206     // SUMMARIZE RESULTS
00207     std::cout << "Number of test queries: " << num_points << std::endl;
00208     std::cout << "Tree build time: " << ( (double)build_time ) / CLOCKS_PER_SEC << " seconds" << std::endl;
00209     std::cout << "Total query time: " << ( (double)test_time ) / CLOCKS_PER_SEC << " seconds" << std::endl;
00210     std::cout << "Time per query: " << ( (double)test_time ) / CLOCKS_PER_SEC / num_points << " seconds" << std::endl;
00211     std::cout << "Tree depth: " << actual_depth << std::endl;
00212     if( actual_depth )
00213         std::cout << "Total query time/depth: " << ( (double)test_time ) / CLOCKS_PER_SEC / actual_depth << " seconds"
00214                   << std::endl;
00215     std::cout << std::endl;
00216     std::cout << "Bytes before tree construction: " << init_total_storage << std::endl;
00217     std::cout << "Bytes after tree construction: " << fini_total_storage << std::endl;
00218     std::cout << "Difference: " << fini_total_storage - init_total_storage << " bytes" << std::endl;
00219     return fail;
00220 }
00221 
00222 void fail( ErrorCode error_code, const char* file, int line )
00223 {
00224     std::cerr << "Internal error (error code " << error_code << ") at " << file << ":" << line << std::endl;
00225     abort();
00226 }
00227 
00228 const int HexSign[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
00229                             { -1, -1, 1 },  { 1, -1, 1 },  { 1, 1, 1 },  { -1, 1, 1 } };
00230 
00231 static CartVect random_point_in_hex( Interface& mb, EntityHandle hex )
00232 {
00233     const double f = RAND_MAX / 2;
00234     CartVect xi( ( (double)rand() - f ) / f, ( (double)rand() - f ) / f, ( (double)rand() - f ) / f );
00235     CartVect coords[8];
00236     const EntityHandle* conn;
00237     int len;
00238     ErrorCode rval = mb.get_connectivity( hex, conn, len, true );
00239     if( len != 8 && MB_SUCCESS != rval )
00240     {
00241         std::cerr << "Invalid element" << std::endl;
00242         assert( false );
00243         abort();
00244     }
00245     rval = mb.get_coords( conn, 8, reinterpret_cast< double* >( coords ) );
00246     CHK( rval );
00247 
00248     CartVect point( 0, 0, 0 );
00249     for( unsigned i = 0; i < 8; ++i )
00250     {
00251         double coeff = 0.125;
00252         for( unsigned j = 0; j < 3; ++j )
00253             coeff *= 1 + HexSign[i][j] * xi[j];
00254         point += coeff * coords[i];
00255     }
00256 
00257     return point;
00258 }
00259 
00260 void generate_random_points( Interface& mb, size_t num_points, std::vector< CartVect >& points,
00261                              std::vector< EntityHandle >& point_elems )
00262 {
00263     Range elems;
00264     ErrorCode rval;
00265     rval = mb.get_entities_by_dimension( 0, 3, elems );
00266     CHK( rval );
00267     if( !elems.all_of_type( MBHEX ) )
00268     {
00269         std::cerr << "Warning: ignoring non-hexahedral elements." << std::endl;
00270         std::pair< Range::iterator, Range::iterator > p = elems.equal_range( MBHEX );
00271         elems.erase( p.second, elems.end() );
00272         elems.erase( elems.begin(), p.first );
00273     }
00274     if( elems.empty() )
00275     {
00276         std::cerr << "Input file contains no hexahedral elements." << std::endl;
00277         exit( 1 );
00278     }
00279 
00280     points.resize( num_points );
00281     point_elems.resize( num_points );
00282     const size_t num_elem = elems.size();
00283     for( size_t i = 0; i < num_points; ++i )
00284     {
00285         size_t offset = 0;
00286         for( size_t x = num_elem; x > 0; x /= RAND_MAX )
00287             offset += rand();
00288         offset %= num_elem;
00289         point_elems[i] = elems[offset];
00290         points[i]      = random_point_in_hex( mb, point_elems[i] );
00291     }
00292 }
00293 
00294 void do_kdtree_test( Interface& mb, int tree_depth, int elem_per_leaf, long num_test,
00295                      const std::vector< CartVect >& points, std::vector< EntityHandle >& point_elems,
00296                      clock_t& build_time, clock_t& test_time, size_t& depth )
00297 {
00298     ErrorCode rval;
00299     clock_t init = clock();
00300     AdaptiveKDTree tool( &mb );
00301     EntityHandle root;
00302     std::ostringstream options;
00303     if( tree_depth > 0 ) options << "MAX_DEPTH=" << tree_depth << ";";
00304     if( elem_per_leaf > 0 ) options << "MAX_PER_LEAF=" << elem_per_leaf << ";";
00305     Range all_hexes;
00306     rval = mb.get_entities_by_type( 0, MBHEX, all_hexes );
00307     CHK( rval );
00308     FileOptions opt( options.str().c_str() );
00309     rval = tool.build_tree( all_hexes, &root, &opt );
00310     CHK( rval );
00311     all_hexes.clear();
00312     build_time = clock() - init;
00313 
00314     EntityHandle leaf;
00315     std::vector< EntityHandle > hexes;
00316     std::vector< EntityHandle >::iterator j;
00317     CartVect coords[8];
00318     const EntityHandle* conn;
00319     int len;
00320     for( long i = 0; i < num_test; ++i )
00321     {
00322         const size_t idx = (size_t)i % points.size();
00323         rval             = tool.point_search( points[idx].array(), leaf, 1.0e-10, 1.0e-6, NULL, &root );
00324         CHK( rval );
00325         hexes.clear();
00326         rval = mb.get_entities_by_handle( leaf, hexes );
00327         CHK( rval );
00328         for( j = hexes.begin(); j != hexes.end(); ++j )
00329         {
00330             rval = mb.get_connectivity( *j, conn, len, true );
00331             CHK( rval );
00332             rval = mb.get_coords( conn, 8, reinterpret_cast< double* >( coords ) );
00333             CHK( rval );
00334             if( GeomUtil::point_in_trilinear_hex( coords, points[idx], 1e-12 ) )
00335             {
00336                 point_elems[idx] = *j;
00337                 break;
00338             }
00339         }
00340         if( j == hexes.end() ) point_elems[idx] = 0;
00341     }
00342 
00343     test_time = clock() - build_time;
00344     double tmp_box[3];
00345     unsigned max_d;
00346     tool.get_info( root, tmp_box, tmp_box, max_d );
00347     depth = max_d;
00348 }
00349 
00350 void do_linear_test( Interface& mb, int, int, long num_test, const std::vector< CartVect >& points,
00351                      std::vector< EntityHandle >& point_elems, clock_t& build_time, clock_t& test_time, size_t& depth )
00352 {
00353     clock_t init = clock();
00354     Range hexes;
00355     ErrorCode rval = mb.get_entities_by_type( 0, MBHEX, hexes );
00356     CHK( rval );
00357     depth = 0;
00358     point_elems.resize( points.size() );
00359     build_time = clock() - init;
00360 
00361     CartVect coords[8];
00362     const EntityHandle* conn;
00363     int len;
00364     for( long i = 0; i < num_test; ++i )
00365     {
00366         const size_t idx = (size_t)i % points.size();
00367         for( Range::iterator h = hexes.begin(); h != hexes.end(); ++h )
00368         {
00369             rval = mb.get_connectivity( *h, conn, len, true );
00370             CHK( rval );
00371             rval = mb.get_coords( conn, 8, reinterpret_cast< double* >( coords ) );
00372             CHK( rval );
00373             if( GeomUtil::point_in_trilinear_hex( coords, points[idx], 1e-12 ) )
00374             {
00375                 point_elems[idx] = *h;
00376                 break;
00377             }
00378         }
00379     }
00380 
00381     test_time = clock() - build_time;
00382 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines