MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 #define IS_BUILDING_MB 00002 #include "moab/Core.hpp" 00003 #include "moab/CartVect.hpp" 00004 #include "moab/OrientedBoxTreeTool.hpp" 00005 #include "moab/OrientedBox.hpp" 00006 #include "Internals.hpp" 00007 #include "moab/Range.hpp" 00008 00009 #include <iostream> 00010 #include <iomanip> 00011 #include <cstdio> 00012 #include <limits> 00013 #include <cstdlib> 00014 #include <ctime> 00015 #if !defined( _MSC_VER ) && !defined( __MINGW32__ ) 00016 #include <unistd.h> 00017 #include <sys/stat.h> 00018 #include <fcntl.h> 00019 #endif 00020 00021 using namespace moab; 00022 00023 std::string clock_to_string( clock_t t ); 00024 std::string mem_to_string( unsigned long mem ); 00025 00026 const int MAX_TAG_VALUE = 20; 00027 const char* const TAG_NAME = "OBB_ID"; 00028 const char* const TREE_TAG = "OBB_ROOT"; 00029 const char* root_tag = TREE_TAG; 00030 00031 ErrorCode get_root( Interface* moab, EntityHandle& root ); 00032 EntityHandle build_tree( Interface* interface, OrientedBoxTreeTool::Settings settings ); 00033 void delete_existing_tree( Interface* interface ); 00034 void print_stats( Interface* interface ); 00035 void tag_triangles( Interface* interface ); 00036 void tag_vertices( Interface* interface ); 00037 void write_tree_blocks( Interface* interface, const char* file ); 00038 00039 static void usage( bool err = true ) 00040 { 00041 std::ostream& s = err ? std::cerr : std::cout; 00042 s << "obb_tree_tool [-s|-S] [-d <int>] [-n <int>] <input file> <output file>" << std::endl 00043 << "obb_tree_tool [-h]" << std::endl; 00044 if( !err ) 00045 { 00046 OrientedBoxTreeTool::Settings st; 00047 s << "Tool to build adaptive kd-Tree from triangles" << std::endl; 00048 s << " -s Use range-based sets for tree nodes" << std::endl 00049 << " -S Use vector-based sets for tree nodes" << std::endl 00050 << " -d <int> Specify maximum depth for tree. Default: " << st.max_depth << std::endl 00051 << " -n <int> Specify maximum entities per leaf. Default: " << st.max_leaf_entities << std::endl 00052 << " -m <real> Specify worst split ratio. Default: " << st.worst_split_ratio << std::endl 00053 << " -M <real> Specify best split ratio. Default: " << st.best_split_ratio << std::endl 00054 << " -t Tag triangles will tree cell number." << std::endl 00055 << " -T Write tree boxes to file." << std::endl 00056 << " -N Specify mesh tag containing tree root. Default: \"" << TREE_TAG << '"' << std::endl 00057 << std::endl; 00058 } 00059 exit( err ); 00060 } 00061 00062 #if defined( _MSC_VER ) || defined( __MINGW32__ ) 00063 static void memory_use( unsigned long long& vsize, unsigned long long& rss ) 00064 { 00065 vsize = rss = 0; 00066 } 00067 #else 00068 static void memory_use( unsigned long long& vsize, unsigned long long& rss ) 00069 { 00070 char buffer[512]; 00071 unsigned long lvsize; 00072 long lrss; 00073 int filp = open( "/proc/self/stat", O_RDONLY ); 00074 ssize_t r = read( filp, buffer, sizeof( buffer ) - 1 ); 00075 close( filp ); 00076 if( r < 0 ) r = 0; 00077 lvsize = lrss = 0; 00078 buffer[r] = '\0'; 00079 sscanf( buffer, 00080 "%*d %*s %*c " // pid command state 00081 "%*d %*d " // ppid pgrp 00082 "%*d %*d %*d " // session tty_nr tpgid 00083 "%*u " // flags 00084 "%*u %*u %*u %*u " // minflt cminflt majflt cmajflt 00085 "%*u %*u %*d %*d " // utime stime cutime cstime 00086 "%*d %*d %*d " // priority nice (unused) 00087 "%*d %*u " // itrealval starttime 00088 "%lu %ld", 00089 &lvsize, &lrss ); 00090 rss = lrss * getpagesize(); 00091 vsize = lvsize; 00092 } 00093 #endif 00094 00095 static int parseint( int& i, int argc, char* argv[] ) 00096 { 00097 char* end; 00098 ++i; 00099 if( i == argc ) 00100 { 00101 std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl; 00102 usage(); 00103 } 00104 00105 int result = strtol( argv[i], &end, 0 ); 00106 if( result < 0 || *end ) 00107 { 00108 std::cerr << "Expected positive integer following '" << argv[i - 1] << "'" << std::endl; 00109 usage(); 00110 } 00111 00112 return result; 00113 } 00114 00115 static double parsedouble( int& i, int argc, char* argv[] ) 00116 { 00117 char* end; 00118 ++i; 00119 if( i == argc ) 00120 { 00121 std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl; 00122 usage(); 00123 } 00124 00125 double result = strtod( argv[i], &end ); 00126 if( result < 0 || *end ) 00127 { 00128 std::cerr << "Expected positive real number following '" << argv[i - 1] << "'" << std::endl; 00129 usage(); 00130 } 00131 00132 return result; 00133 } 00134 00135 int main( int argc, char* argv[] ) 00136 { 00137 const char* input_file = 0; 00138 const char* output_file = 0; 00139 const char* tree_file = 0; 00140 OrientedBoxTreeTool::Settings settings; 00141 bool tag_tris = false; 00142 clock_t load_time, build_time, stat_time, tag_time, write_time, block_time; 00143 00144 for( int i = 1; i < argc; ++i ) 00145 { 00146 if( argv[i][0] != '-' ) 00147 { 00148 if( !input_file ) 00149 input_file = argv[i]; 00150 else if( !output_file ) 00151 output_file = argv[i]; 00152 else 00153 usage(); 00154 continue; 00155 } 00156 00157 if( !argv[i][1] || argv[i][2] ) usage(); 00158 00159 switch( argv[i][1] ) 00160 { 00161 case 's': 00162 settings.set_options = MESHSET_SET; 00163 break; 00164 case 'S': 00165 settings.set_options = MESHSET_ORDERED; 00166 break; 00167 case 'd': 00168 settings.max_depth = parseint( i, argc, argv ); 00169 break; 00170 case 'n': 00171 settings.max_leaf_entities = parseint( i, argc, argv ); 00172 break; 00173 case 'm': 00174 settings.worst_split_ratio = parsedouble( i, argc, argv ); 00175 break; 00176 case 'M': 00177 settings.best_split_ratio = parsedouble( i, argc, argv ); 00178 break; 00179 case 't': 00180 tag_tris = true; 00181 break; 00182 case 'T': 00183 if( ++i == argc ) usage(); 00184 tree_file = argv[i]; 00185 break; 00186 case 'N': 00187 if( ++i == argc ) usage(); 00188 root_tag = argv[i]; 00189 break; 00190 case 'h': 00191 usage( false ); 00192 default: 00193 usage(); 00194 } 00195 } 00196 00197 if( !output_file ) usage(); 00198 00199 ErrorCode rval; 00200 Core moab_core; 00201 Interface* interface = &moab_core; 00202 00203 load_time = clock(); 00204 rval = interface->load_mesh( input_file ); 00205 if( MB_SUCCESS != rval ) 00206 { 00207 std::cerr << "Error reading file: " << input_file << std::endl; 00208 exit( 2 ); 00209 } 00210 load_time = clock() - load_time; 00211 00212 delete_existing_tree( interface ); 00213 00214 std::cout << "Building tree..." << std::endl; 00215 build_time = clock(); 00216 build_tree( interface, settings ); 00217 build_time = clock() - build_time; 00218 00219 std::cout << "Calculating stats..." << std::endl; 00220 print_stats( interface ); 00221 stat_time = clock() - build_time; 00222 00223 if( tag_tris ) 00224 { 00225 std::cout << "Tagging tree..." << std::endl; 00226 tag_triangles( interface ); 00227 tag_vertices( interface ); 00228 } 00229 tag_time = clock() - stat_time; 00230 00231 std::cout << "Writing file... "; 00232 std::cout.flush(); 00233 rval = interface->write_mesh( output_file ); 00234 if( MB_SUCCESS != rval ) 00235 { 00236 std::cerr << "Error writing file: " << output_file << std::endl; 00237 exit( 3 ); 00238 } 00239 write_time = clock() - tag_time; 00240 std::cout << "Wrote " << output_file << std::endl; 00241 00242 if( tree_file ) 00243 { 00244 std::cout << "Writing tree block rep..."; 00245 std::cout.flush(); 00246 write_tree_blocks( interface, tree_file ); 00247 std::cout << "Wrote " << tree_file << std::endl; 00248 } 00249 block_time = clock() - write_time; 00250 00251 std::cout << "Times: " 00252 << " Load" 00253 << " Build" 00254 << " Stats" 00255 << " Write"; 00256 if( tag_tris ) std::cout << "Tag Sets"; 00257 if( tree_file ) std::cout << "Block "; 00258 std::cout << std::endl; 00259 00260 std::cout << " " << std::setw( 8 ) << clock_to_string( load_time ) << std::setw( 8 ) 00261 << clock_to_string( build_time ) << std::setw( 8 ) << clock_to_string( stat_time ) << std::setw( 8 ) 00262 << clock_to_string( write_time ); 00263 if( tag_tris ) std::cout << std::setw( 8 ) << clock_to_string( tag_time ); 00264 if( tree_file ) std::cout << std::setw( 8 ) << clock_to_string( block_time ); 00265 std::cout << std::endl; 00266 00267 return 0; 00268 } 00269 00270 ErrorCode get_root( Interface* moab, EntityHandle& root ) 00271 { 00272 Tag tag; 00273 ErrorCode rval; 00274 00275 rval = moab->tag_get_handle( root_tag, 1, MB_TYPE_HANDLE, tag ); 00276 if( MB_SUCCESS != rval ) return rval; 00277 00278 const EntityHandle mesh = 0; 00279 return moab->tag_get_data( tag, &mesh, 1, &root ); 00280 } 00281 00282 void delete_existing_tree( Interface* interface ) 00283 { 00284 EntityHandle root; 00285 ErrorCode rval = get_root( interface, root ); 00286 if( MB_SUCCESS == rval ) 00287 { 00288 OrientedBoxTreeTool tool( interface ); 00289 rval = tool.delete_tree( root ); 00290 if( MB_SUCCESS != rval ) 00291 { 00292 std::cerr << "Failed to destroy existing trees. Aborting" << std::endl; 00293 exit( 5 ); 00294 } 00295 } 00296 } 00297 00298 EntityHandle build_tree( Interface* interface, OrientedBoxTreeTool::Settings settings ) 00299 { 00300 ErrorCode rval; 00301 EntityHandle root = 0; 00302 Range triangles; 00303 00304 rval = interface->get_entities_by_type( 0, MBTRI, triangles ); 00305 if( MB_SUCCESS != rval || triangles.empty() ) 00306 { 00307 std::cerr << "No triangles from which to build tree." << std::endl; 00308 exit( 4 ); 00309 } 00310 00311 OrientedBoxTreeTool tool( interface ); 00312 rval = tool.build( triangles, root, &settings ); 00313 if( MB_SUCCESS != rval || !root ) 00314 { 00315 std::cerr << "Tree construction failed." << std::endl; 00316 exit( 4 ); 00317 } 00318 00319 // store tree root 00320 Tag roottag; 00321 rval = interface->tag_get_handle( root_tag, 1, MB_TYPE_HANDLE, roottag, MB_TAG_CREAT | MB_TAG_SPARSE ); 00322 if( MB_SUCCESS != rval ) 00323 { 00324 std::cout << "Failed to create root tag: \"" << root_tag << '"' << std::endl; 00325 exit( 2 ); 00326 } 00327 const EntityHandle mesh = 0; 00328 rval = interface->tag_set_data( roottag, &mesh, 1, &root ); 00329 if( MB_SUCCESS != rval ) 00330 { 00331 std::cout << "Failed to set root tag: \"" << root_tag << '"' << std::endl; 00332 exit( 2 ); 00333 } 00334 00335 return root; 00336 } 00337 00338 std::string clock_to_string( clock_t t ) 00339 { 00340 char unit[5] = "s"; 00341 char buffer[256]; 00342 double dt = t / (double)CLOCKS_PER_SEC; 00343 // if (dt > 300) { 00344 // dt /= 60; 00345 // strcpy( unit, "min" ); 00346 //} 00347 // if (dt > 300) { 00348 // dt /= 60; 00349 // strcpy( unit, "hr" ); 00350 //} 00351 // if (dt > 100) { 00352 // dt /= 24; 00353 // strcpy( unit, "days" ); 00354 //} 00355 sprintf( buffer, "%0.2f%s", dt, unit ); 00356 return buffer; 00357 } 00358 00359 std::string mem_to_string( unsigned long mem ) 00360 { 00361 char unit[3] = "B"; 00362 if( mem > 9 * 1024 ) 00363 { 00364 mem = ( mem + 512 ) / 1024; 00365 strcpy( unit, "kB" ); 00366 } 00367 if( mem > 9 * 1024 ) 00368 { 00369 mem = ( mem + 512 ) / 1024; 00370 strcpy( unit, "MB" ); 00371 } 00372 if( mem > 9 * 1024 ) 00373 { 00374 mem = ( mem + 512 ) / 1024; 00375 strcpy( unit, "GB" ); 00376 } 00377 char buffer[256]; 00378 sprintf( buffer, "%lu %s", mem, unit ); 00379 return buffer; 00380 } 00381 00382 template < typename T > 00383 struct SimpleStat 00384 { 00385 T min, max, sum, sqr; 00386 size_t count; 00387 SimpleStat(); 00388 void add( T value ); 00389 double avg() const 00390 { 00391 return (double)sum / count; 00392 } 00393 double rms() const 00394 { 00395 return sqrt( (double)sqr / count ); 00396 } 00397 double dev() const 00398 { 00399 return sqrt( ( count * (double)sqr - (double)sum * (double)sum ) / ( (double)count * ( count - 1 ) ) ); 00400 } 00401 }; 00402 00403 template < typename T > 00404 SimpleStat< T >::SimpleStat() 00405 : min( std::numeric_limits< T >::max() ), max( std::numeric_limits< T >::min() ), sum( 0 ), sqr( 0 ), count( 0 ) 00406 { 00407 } 00408 00409 void print_stats( Interface* interface ) 00410 { 00411 EntityHandle root; 00412 Range range; 00413 ErrorCode rval; 00414 rval = get_root( interface, root ); 00415 if( MB_SUCCESS != rval ) 00416 { 00417 std::cerr << "Internal error: Failed to retrieve root." << std::endl; 00418 exit( 5 ); 00419 } 00420 OrientedBoxTreeTool tool( interface ); 00421 00422 Range tree_sets, triangles, verts; 00423 // interface->get_child_meshsets( root, tree_sets, 0 ); 00424 interface->get_entities_by_type( 0, MBENTITYSET, tree_sets ); 00425 tree_sets.erase( tree_sets.begin(), Range::lower_bound( tree_sets.begin(), tree_sets.end(), root ) ); 00426 interface->get_entities_by_type( 0, MBTRI, triangles ); 00427 interface->get_entities_by_type( 0, MBVERTEX, verts ); 00428 triangles.merge( verts ); 00429 tree_sets.insert( root ); 00430 unsigned long long set_used, set_amortized, set_store_used, set_store_amortized, set_tag_used, set_tag_amortized, 00431 tri_used, tri_amortized; 00432 interface->estimated_memory_use( tree_sets, &set_used, &set_amortized, &set_store_used, &set_store_amortized, 0, 0, 00433 0, 0, &set_tag_used, &set_tag_amortized ); 00434 interface->estimated_memory_use( triangles, &tri_used, &tri_amortized ); 00435 00436 int num_tri = 0; 00437 interface->get_number_entities_by_type( 0, MBTRI, num_tri ); 00438 00439 tool.stats( root, std::cout ); 00440 00441 unsigned long long real_rss, real_vsize; 00442 memory_use( real_vsize, real_rss ); 00443 00444 printf( "------------------------------------------------------------------\n" ); 00445 printf( "\nmemory: used amortized\n" ); 00446 printf( " ---------- ----------\n" ); 00447 printf( "triangles %10s %10s\n", mem_to_string( tri_used ).c_str(), mem_to_string( tri_amortized ).c_str() ); 00448 printf( "sets (total)%10s %10s\n", mem_to_string( set_used ).c_str(), mem_to_string( set_amortized ).c_str() ); 00449 printf( "sets %10s %10s\n", mem_to_string( set_store_used ).c_str(), 00450 mem_to_string( set_store_amortized ).c_str() ); 00451 printf( "set tags %10s %10s\n", mem_to_string( set_tag_used ).c_str(), 00452 mem_to_string( set_tag_amortized ).c_str() ); 00453 printf( "total real %10s %10s\n", mem_to_string( real_rss ).c_str(), mem_to_string( real_vsize ).c_str() ); 00454 printf( "------------------------------------------------------------------\n" ); 00455 } 00456 00457 static int hash_handle( EntityHandle handle ) 00458 { 00459 EntityID h = ID_FROM_HANDLE( handle ); 00460 return (int)( ( h * 13 + 7 ) % MAX_TAG_VALUE ) + 1; 00461 } 00462 00463 class TriTagger : public OrientedBoxTreeTool::Op 00464 { 00465 private: 00466 Interface* mMB; 00467 Tag mTag; 00468 std::vector< EntityHandle > mHandles; 00469 std::vector< int > mTagData; 00470 00471 public: 00472 TriTagger( Tag tag, Interface* moab ) : mMB( moab ), mTag( tag ) {} 00473 00474 ErrorCode visit( EntityHandle, int, bool& descent ) 00475 { 00476 descent = true; 00477 return MB_SUCCESS; 00478 } 00479 00480 ErrorCode leaf( EntityHandle node ) 00481 { 00482 mHandles.clear(); 00483 mMB->get_entities_by_handle( node, mHandles ); 00484 mTagData.clear(); 00485 mTagData.resize( mHandles.size(), hash_handle( node ) ); 00486 mMB->tag_set_data( mTag, &mHandles[0], mHandles.size(), &mTagData[0] ); 00487 return MB_SUCCESS; 00488 } 00489 }; 00490 00491 void tag_triangles( Interface* moab ) 00492 { 00493 EntityHandle root; 00494 ErrorCode rval = get_root( moab, root ); 00495 if( MB_SUCCESS != rval ) 00496 { 00497 std::cerr << "Internal error: Failed to retrieve tree." << std::endl; 00498 exit( 5 ); 00499 } 00500 00501 Tag tag; 00502 int zero = 0; 00503 moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero ); 00504 TriTagger op( tag, moab ); 00505 00506 OrientedBoxTreeTool tool( moab ); 00507 rval = tool.preorder_traverse( root, op ); 00508 if( MB_SUCCESS != rval ) 00509 { 00510 std::cerr << "Internal error tagging triangles" << std::endl; 00511 exit( 5 ); 00512 } 00513 } 00514 00515 class VtxTagger : public OrientedBoxTreeTool::Op 00516 { 00517 private: 00518 Interface* mMB; 00519 Tag mTag; 00520 std::vector< EntityHandle > mHandles; 00521 std::vector< EntityHandle > mConn; 00522 std::vector< int > mTagData; 00523 00524 public: 00525 VtxTagger( Tag tag, Interface* moab ) : mMB( moab ), mTag( tag ) {} 00526 00527 ErrorCode visit( EntityHandle, int, bool& descent ) 00528 { 00529 descent = true; 00530 return MB_SUCCESS; 00531 } 00532 00533 ErrorCode leaf( EntityHandle node ) 00534 { 00535 mHandles.clear(); 00536 mMB->get_entities_by_handle( node, mHandles ); 00537 mConn.clear(); 00538 mMB->get_connectivity( &mHandles[0], mHandles.size(), mConn ); 00539 mTagData.clear(); 00540 mTagData.resize( mConn.size(), hash_handle( node ) ); 00541 mMB->tag_set_data( mTag, &mConn[0], mConn.size(), &mTagData[0] ); 00542 return MB_SUCCESS; 00543 } 00544 }; 00545 00546 void tag_vertices( Interface* moab ) 00547 { 00548 EntityHandle root; 00549 ErrorCode rval = get_root( moab, root ); 00550 if( MB_SUCCESS != rval ) 00551 { 00552 std::cerr << "Internal error: Failed to retrieve tree." << std::endl; 00553 exit( 5 ); 00554 } 00555 00556 Tag tag; 00557 int zero = 0; 00558 moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero ); 00559 VtxTagger op( tag, moab ); 00560 00561 OrientedBoxTreeTool tool( moab ); 00562 rval = tool.preorder_traverse( root, op ); 00563 if( MB_SUCCESS != rval ) 00564 { 00565 std::cerr << "Internal error tagging vertices" << std::endl; 00566 exit( 5 ); 00567 } 00568 } 00569 00570 class LeafHexer : public OrientedBoxTreeTool::Op 00571 { 00572 private: 00573 OrientedBoxTreeTool* mTool; 00574 Interface* mOut; 00575 Tag mTag; 00576 std::vector< EntityHandle > mHandles; 00577 std::vector< EntityHandle > mConn; 00578 std::vector< int > mTagData; 00579 00580 public: 00581 LeafHexer( OrientedBoxTreeTool* tool, Interface* mb2, Tag tag ) : mTool( tool ), mOut( mb2 ), mTag( tag ) {} 00582 00583 ErrorCode visit( EntityHandle, int, bool& descent ) 00584 { 00585 descent = true; 00586 return MB_SUCCESS; 00587 } 00588 00589 ErrorCode leaf( EntityHandle node ) 00590 { 00591 OrientedBox box; 00592 ErrorCode rval = mTool->box( node, box ); 00593 if( MB_SUCCESS != rval ) return rval; 00594 EntityHandle h; 00595 rval = box.make_hex( h, mOut ); 00596 if( MB_SUCCESS != rval ) return rval; 00597 int i = hash_handle( node ); 00598 return mOut->tag_set_data( mTag, &h, 1, &i ); 00599 } 00600 }; 00601 00602 void write_tree_blocks( Interface* interface, const char* file ) 00603 { 00604 EntityHandle root; 00605 ErrorCode rval = get_root( interface, root ); 00606 if( MB_SUCCESS != rval ) 00607 { 00608 std::cerr << "Internal error: Failed to retrieve tree." << std::endl; 00609 exit( 5 ); 00610 } 00611 00612 Core moab2; 00613 Tag tag; 00614 int zero = 0; 00615 moab2.tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero ); 00616 00617 OrientedBoxTreeTool tool( interface ); 00618 LeafHexer op( &tool, &moab2, tag ); 00619 rval = tool.preorder_traverse( root, op ); 00620 if( MB_SUCCESS != rval ) 00621 { 00622 std::cerr << "Internal error: failed to construct leaf hexes" << std::endl; 00623 exit( 5 ); 00624 } 00625 00626 moab2.write_mesh( file ); 00627 }