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