MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 #include <iostream> 00002 #include <vector> 00003 #include <string> 00004 #include <iomanip> 00005 #include <fstream> 00006 #if defined( __MINGW32__ ) 00007 #include <sys/time.h> 00008 #else 00009 #include <ctime> 00010 #endif 00011 00012 #include <ctime> 00013 #include <cmath> 00014 #include <cassert> 00015 #include <cfloat> 00016 00017 #include "moab/Core.hpp" 00018 #include "moab/Interface.hpp" 00019 #include "moab/verdict/VerdictWrapper.hpp" 00020 #include "moab/NestedRefine.hpp" 00021 00022 #ifdef MOAB_HAVE_MPI 00023 #include "moab_mpi.h" 00024 #include "moab/ParallelComm.hpp" 00025 #include "MBParallelConventions.h" 00026 #endif 00027 00028 /* Exit values */ 00029 #define SUCCESS 0 00030 #define USAGE_ERROR 1 00031 #define NOT_IMPLEMENTED 2 00032 00033 using namespace moab; 00034 00035 static void print_usage( const char* name, std::ostream& stream ) 00036 { 00037 stream << "Usage: " << name << " <options> <input_file> [<input_file2> ...]" << std::endl 00038 << "Options: " << std::endl 00039 << "\t-h - Print this help text and exit." << std::endl 00040 << "\t-d - Dimension of the mesh." << std::endl 00041 << "\t-n - Exact or a maximum number of levels for the hierarchy. Default 1." << std::endl 00042 << "\t-L - Degree of refinement for each level. Pass an array or a number. It " 00043 "is mandatory to pass dimension and num_levels before to use this option. If this flag " 00044 "is not used, a default of deg 2 refinement is used. " 00045 << std::endl 00046 << "\t-V - Pass a desired volume (absolute) . This will generate a hierarchy " 00047 "such that the maximum volume is reduced to the given amount approximately. The length " 00048 "of the hierarchy can be constrained further if a maximum number of levels is passed. " 00049 "It is mandatory to pass the dimension for this option. " 00050 << std::endl 00051 << "\t-q - Prints out the maximum volume of the mesh and exits the program. " 00052 "This option can be used as a guide to volume constrainted mesh hierarchy later. " 00053 << std::endl 00054 << "\t-t - Print out the time taken to generate hierarchy." << std::endl 00055 << "\t-s - Print out the mesh sizes of each level of the generated hierarchy." << std::endl 00056 << "\t-o - Specify true for output files for the mesh levels of the hierarchy." << std::endl 00057 //<< "\t-O option - Specify read option." << std::endl 00058 #ifdef MOAB_HAVE_MPI 00059 << "\t-p[0|1|2] - Read in parallel[0], optionally also doing resolve_shared_ents (1) " 00060 "and exchange_ghosts (2)" 00061 << std::endl 00062 #endif 00063 ; 00064 exit( USAGE_ERROR ); 00065 } 00066 00067 static void usage_error( const char* name ) 00068 { 00069 print_usage( name, std::cerr ); 00070 #ifdef MOAB_HAVE_MPI 00071 MPI_Finalize(); 00072 #endif 00073 exit( USAGE_ERROR ); 00074 } 00075 00076 bool parse_id_list( const char* string, int dim, int nval, std::vector< int >& results ); 00077 00078 bool make_opts_string( std::vector< std::string > options, std::string& opts ); 00079 00080 ErrorCode get_degree_seq( Core& mb, EntityHandle fileset, int dim, double desired_vol, int& num_levels, 00081 std::vector< int >& level_degs ); 00082 00083 ErrorCode get_max_volume( Core& mb, EntityHandle fileset, int dim, double& vmax ); 00084 00085 int main( int argc, char* argv[] ) 00086 { 00087 int proc_id = 0, size = 1; 00088 #ifdef MOAB_HAVE_MPI 00089 MPI_Init( &argc, &argv ); 00090 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id ); 00091 MPI_Comm_size( MPI_COMM_WORLD, &size ); 00092 #endif 00093 00094 int num_levels = 0, dim = 0; 00095 std::vector< int > level_degrees; 00096 bool optimize = false; 00097 bool do_flag = true; 00098 bool print_times = false, print_size = false, output = false; 00099 bool parallel = false, resolve_shared = false, exchange_ghosts = false; 00100 bool printusage = false, parselevels = true; 00101 bool qc_vol = false, only_quality = false; 00102 double cvol = 0; 00103 std::string infile; 00104 00105 int i = 1; 00106 while( i < argc ) 00107 { 00108 if( !argv[i][0] && 0 == proc_id ) 00109 { 00110 usage_error( argv[0] ); 00111 #ifdef MOAB_HAVE_MPI 00112 MPI_Finalize(); 00113 #endif 00114 } 00115 00116 if( do_flag && argv[i][0] == '-' ) 00117 { 00118 switch( argv[i][1] ) 00119 { 00120 // do flag arguments: 00121 case '-': 00122 do_flag = false; 00123 break; 00124 case 't': 00125 print_times = true; 00126 break; 00127 case 's': 00128 print_size = true; 00129 break; 00130 case 'h': 00131 case 'H': 00132 print_usage( argv[0], std::cerr ); 00133 printusage = true; 00134 break; 00135 case 'd': 00136 dim = atol( argv[i + 1] ); 00137 ++i; 00138 break; 00139 case 'n': 00140 num_levels = atol( argv[i + 1] ); 00141 ++i; 00142 break; 00143 case 'L': 00144 if( dim != 0 && num_levels != 0 ) 00145 { 00146 parselevels = parse_id_list( argv[i + 1], dim, num_levels, level_degrees ); 00147 ++i; 00148 } 00149 else 00150 { 00151 print_usage( argv[0], std::cerr ); 00152 printusage = true; 00153 } 00154 break; 00155 case 'V': 00156 qc_vol = true; 00157 cvol = strtod( argv[i + 1], NULL ); 00158 ++i; 00159 break; 00160 case 'q': 00161 only_quality = true; 00162 break; 00163 case 'o': 00164 output = true; 00165 break; 00166 case 'O': 00167 optimize = true; 00168 break; 00169 #ifdef MOAB_HAVE_MPI 00170 case 'p': 00171 parallel = true; 00172 resolve_shared = true; 00173 if( argv[i][1] == '1' ) exchange_ghosts = true; 00174 break; 00175 #endif 00176 default: 00177 ++i; 00178 switch( argv[i - 1][1] ) 00179 { 00180 // case 'O': read_opts.push_back(argv[i]); break; 00181 default: 00182 std::cerr << "Invalid option: " << argv[i] << std::endl; 00183 } 00184 } 00185 ++i; 00186 } 00187 // do file names 00188 else 00189 { 00190 infile = argv[i]; 00191 ++i; 00192 } 00193 } 00194 00195 if( infile.empty() && !printusage ) print_usage( argv[0], std::cerr ); 00196 00197 if( !parselevels ) exit( USAGE_ERROR ); 00198 00199 #ifdef MOAB_HAVE_MPI 00200 parallel = true; 00201 resolve_shared = true; 00202 #endif 00203 00204 ErrorCode error; 00205 00206 Core* moab = new Core; 00207 Interface* mb = (Interface*)moab; 00208 EntityHandle fileset; 00209 00210 // Create a fileset 00211 error = mb->create_meshset( MESHSET_SET, fileset );MB_CHK_ERR( error ); 00212 00213 // Set the read options for parallel file loading 00214 std::vector< std::string > read_opts, write_opts; 00215 std::string read_options, write_options; 00216 00217 if( parallel && size > 1 ) 00218 { 00219 read_opts.push_back( "PARALLEL=READ_PART" ); 00220 read_opts.push_back( "PARTITION=PARALLEL_PARTITION" ); 00221 if( resolve_shared ) read_opts.push_back( "PARALLEL_RESOLVE_SHARED_ENTS" ); 00222 if( exchange_ghosts ) read_opts.push_back( "PARALLEL_GHOSTS=3.0.1" ); 00223 00224 /* if (output) 00225 { 00226 write_opts.push_back(";;PARALLEL=WRITE_PART"); 00227 std::cout<<"Here"<<std::endl; 00228 }*/ 00229 } 00230 00231 if( !make_opts_string( read_opts, read_options ) || !make_opts_string( write_opts, write_options ) ) 00232 { 00233 #ifdef MOAB_HAVE_MPI 00234 MPI_Finalize(); 00235 #endif 00236 return USAGE_ERROR; 00237 } 00238 00239 // Load file 00240 std::cout << "READ OPTS=" << (char*)read_options.c_str() << std::endl; 00241 error = mb->load_file( infile.c_str(), &fileset, read_options.c_str() );MB_CHK_ERR( error ); 00242 00243 // Create the nestedrefine instance 00244 00245 #ifdef MOAB_HAVE_MPI 00246 ParallelComm* pc = new ParallelComm( mb, MPI_COMM_WORLD ); 00247 NestedRefine* uref = new NestedRefine( moab, pc, fileset ); 00248 #else 00249 NestedRefine* uref = new NestedRefine( moab ); 00250 #endif 00251 00252 std::vector< EntityHandle > lsets; 00253 00254 // std::cout<<"Read input file"<<std::endl; 00255 00256 if( only_quality ) 00257 { 00258 double vmax; 00259 error = get_max_volume( *moab, fileset, dim, vmax );MB_CHK_ERR( error ); 00260 #ifdef MOAB_HAVE_MPI 00261 int rank = 0; 00262 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00263 if( rank == 0 ) std::cout << "Maximum global volume = " << vmax << std::endl; 00264 #else 00265 std::cout << "Maximum volume = " << vmax << std::endl; 00266 #endif 00267 exit( SUCCESS ); 00268 } 00269 00270 // If a volume constraint is given, find an optimal degree sequence to reach the desired volume 00271 // constraint. 00272 if( qc_vol ) 00273 { 00274 error = get_degree_seq( *moab, fileset, dim, cvol, num_levels, level_degrees );MB_CHK_ERR( error ); 00275 00276 if( dim == 0 ) print_usage( argv[0], std::cerr ); 00277 } 00278 00279 if( num_levels == 0 ) num_levels = 1; 00280 00281 int* ldeg = new int[num_levels]; 00282 // std::cout<<"level_degrees.size = "<<level_degrees.size()<<std::endl; 00283 if( level_degrees.empty() ) 00284 { 00285 for( int l = 0; l < num_levels; l++ ) 00286 ldeg[l] = 2; 00287 } 00288 else 00289 { 00290 for( int l = 0; l < num_levels; l++ ) 00291 ldeg[l] = level_degrees[l]; 00292 } 00293 00294 std::cout << "Starting hierarchy generation" << std::endl; 00295 00296 std::cout << "opt = " << optimize << std::endl; 00297 00298 error = uref->generate_mesh_hierarchy( num_levels, ldeg, lsets, optimize );MB_CHK_ERR( error ); 00299 00300 if( print_times ) 00301 { 00302 std::cout << "Finished hierarchy generation in " << uref->timeall.tm_total << " secs" << std::endl; 00303 if( parallel ) 00304 { 00305 std::cout << "Time taken for refinement " << uref->timeall.tm_refine << " secs" << std::endl; 00306 std::cout << "Time taken for resolving shared interface " << uref->timeall.tm_resolve << " secs" 00307 << std::endl; 00308 } 00309 } 00310 else 00311 std::cout << "Finished hierarchy generation " << std::endl; 00312 00313 if( print_size ) 00314 { 00315 Range all_ents, ents[4]; 00316 error = mb->get_entities_by_handle( fileset, all_ents );MB_CHK_ERR( error ); 00317 00318 for( int k = 0; k < 4; k++ ) 00319 ents[k] = all_ents.subset_by_dimension( k ); 00320 00321 std::cout << std::endl; 00322 if( qc_vol ) 00323 { 00324 double volume; 00325 error = get_max_volume( *moab, fileset, dim, volume );MB_CHK_ERR( error ); 00326 std::cout << "Mesh size for level 0" 00327 << " :: nverts = " << ents[0].size() << ", nedges = " << ents[1].size() 00328 << ", nfaces = " << ents[2].size() << ", ncells = " << ents[3].size() << " :: Vmax = " << volume 00329 << std::endl; 00330 } 00331 else 00332 std::cout << "Mesh size for level 0" 00333 << " :: nverts = " << ents[0].size() << ", nedges = " << ents[1].size() 00334 << ", nfaces = " << ents[2].size() << ", ncells = " << ents[3].size() << std::endl; 00335 00336 for( int l = 0; l < num_levels; l++ ) 00337 { 00338 all_ents.clear(); 00339 ents[0].clear(); 00340 ents[1].clear(); 00341 ents[2].clear(); 00342 ents[3].clear(); 00343 error = mb->get_entities_by_handle( lsets[l + 1], all_ents );MB_CHK_ERR( error ); 00344 00345 for( int k = 0; k < 4; k++ ) 00346 ents[k] = all_ents.subset_by_dimension( k ); 00347 00348 // std::cout<<std::endl; 00349 00350 if( qc_vol ) 00351 { 00352 double volume; 00353 error = get_max_volume( *moab, lsets[l + 1], dim, volume );MB_CHK_ERR( error ); 00354 std::cout << "Mesh size for level " << l + 1 << " :: nverts = " << ents[0].size() 00355 << ", nedges = " << ents[1].size() << ", nfaces = " << ents[2].size() 00356 << ", ncells = " << ents[3].size() << " :: Vmax = " << volume << std::endl; 00357 } 00358 else 00359 std::cout << "Mesh size for level " << l + 1 << " :: nverts = " << ents[0].size() 00360 << ", nedges = " << ents[1].size() << ", nfaces = " << ents[2].size() 00361 << ", ncells = " << ents[3].size() << std::endl; 00362 } 00363 } 00364 00365 if( output ) 00366 { 00367 for( int l = 0; l < num_levels; l++ ) 00368 { 00369 std::string::size_type idx1 = infile.find_last_of( "\\/" ); 00370 std::string::size_type idx2 = infile.find_last_of( "." ); 00371 std::string file = infile.substr( idx1 + 1, idx2 - idx1 - 1 ); 00372 std::stringstream out; 00373 if( parallel ) 00374 out << "_ML_" << l + 1 << ".h5m"; 00375 else 00376 out << "_ML_" << l + 1 << ".vtk"; 00377 file = file + out.str(); 00378 const char* output_file = file.c_str(); 00379 #ifdef MOAB_HAVE_MPI 00380 error = mb->write_file( output_file, 0, ";;PARALLEL=WRITE_PART", &lsets[l + 1], 1 );MB_CHK_ERR( error ); 00381 #else 00382 error = mb->write_file( output_file, 0, NULL, &lsets[l + 1], 1 );MB_CHK_ERR( error ); 00383 #endif 00384 // const char* output_file = file.c_str(); 00385 // error = mb->write_file(output_file, 0, write_options.c_str(), &lsets[l+1], 00386 // 1);MB_CHK_ERR(error); 00387 // mb->list_entity(lsets[l+1]); 00388 // mb->write_file(output_file, 0, "PARALLEL=WRITE_PART", &lsets[l+1], 1); 00389 } 00390 } 00391 00392 delete uref; 00393 #ifdef MOAB_HAVE_MPI 00394 delete pc; 00395 #endif 00396 00397 delete[] ldeg; 00398 delete moab; 00399 00400 #ifdef MOAB_HAVE_MPI 00401 MPI_Finalize(); 00402 #endif 00403 00404 exit( SUCCESS ); 00405 } 00406 00407 ErrorCode get_degree_seq( Core& mb, EntityHandle fileset, int dim, double desired_vol, int& num_levels, 00408 std::vector< int >& level_degs ) 00409 { 00410 // Find max volume 00411 double vmax_global; 00412 ErrorCode error = get_max_volume( mb, fileset, dim, vmax_global );MB_CHK_ERR( error ); 00413 00414 int init_nl = num_levels; 00415 num_levels = 0; 00416 level_degs.clear(); 00417 00418 // Find a sequence that reduces the maximum volume to desired. 00419 // desired_vol should be a fraction or absolute value ? 00420 00421 // double remV = vmax_global*desired_vol/dim; 00422 double Vdesired = desired_vol; 00423 double try_x; 00424 double remV = vmax_global; 00425 int degs[3][3] = { { 5, 3, 2 }, { 25, 9, 4 }, { 0, 27, 8 } }; 00426 00427 if( dim == 1 || dim == 2 ) 00428 { 00429 while( remV - Vdesired >= 0 ) 00430 { 00431 try_x = degs[dim - 1][0]; 00432 if( ( remV / try_x - Vdesired ) >= 0 ) 00433 { 00434 level_degs.push_back( 5 ); 00435 num_levels += 1; 00436 remV /= try_x; 00437 } 00438 else 00439 { 00440 try_x = degs[dim - 1][1]; 00441 if( ( remV / try_x - Vdesired ) >= 0 ) 00442 { 00443 level_degs.push_back( 3 ); 00444 num_levels += 1; 00445 remV /= try_x; 00446 } 00447 else 00448 { 00449 try_x = degs[dim - 1][2]; 00450 if( ( remV / try_x - Vdesired ) >= 0 ) 00451 { 00452 level_degs.push_back( 2 ); 00453 num_levels += 1; 00454 remV /= try_x; 00455 } 00456 else 00457 break; 00458 } 00459 } 00460 } 00461 } 00462 else 00463 { 00464 while( remV - Vdesired >= 0 ) 00465 { 00466 try_x = degs[dim - 1][1]; 00467 if( ( remV / try_x - Vdesired ) >= 0 ) 00468 { 00469 level_degs.push_back( 3 ); 00470 num_levels += 1; 00471 remV /= try_x; 00472 } 00473 else 00474 { 00475 try_x = degs[dim - 1][2]; 00476 if( ( remV / try_x - Vdesired ) >= 0 ) 00477 { 00478 level_degs.push_back( 2 ); 00479 num_levels += 1; 00480 remV /= try_x; 00481 } 00482 else 00483 break; 00484 } 00485 } 00486 } 00487 00488 if( init_nl != 0 && init_nl < num_levels ) 00489 { 00490 for( int i = level_degs.size(); i >= init_nl; i-- ) 00491 level_degs.pop_back(); 00492 num_levels = init_nl; 00493 } 00494 00495 return MB_SUCCESS; 00496 } 00497 00498 ErrorCode get_max_volume( Core& mb, EntityHandle fileset, int dim, double& vmax ) 00499 { 00500 ErrorCode error; 00501 VerdictWrapper vw( &mb ); 00502 QualityType q; 00503 00504 switch( dim ) 00505 { 00506 case 1: 00507 q = MB_LENGTH; 00508 break; 00509 case 2: 00510 q = MB_AREA; 00511 break; 00512 case 3: 00513 q = MB_VOLUME; 00514 break; 00515 default: 00516 return MB_FAILURE; 00517 break; 00518 } 00519 00520 // Get all entities of the highest dimension which is passed as a command line argument. 00521 Range allents, owned; 00522 error = mb.get_entities_by_handle( fileset, allents );MB_CHK_ERR( error ); 00523 owned = allents.subset_by_dimension( dim );MB_CHK_ERR( error ); 00524 00525 // Get all owned entities 00526 #ifdef MOAB_HAVE_MPI 00527 int size = 1; 00528 MPI_Comm_size( MPI_COMM_WORLD, &size ); 00529 int mpi_err; 00530 if( size > 1 ) 00531 { 00532 // filter the entities not owned, so we do not process them more than once 00533 ParallelComm* pcomm = moab::ParallelComm::get_pcomm( &mb, 0 ); 00534 Range current = owned; 00535 owned.clear(); 00536 error = pcomm->filter_pstatus( current, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned ); 00537 if( error != MB_SUCCESS ) 00538 { 00539 MPI_Finalize(); 00540 return MB_FAILURE; 00541 } 00542 } 00543 #endif 00544 00545 double vmax_local = 0; 00546 // Find the maximum volume of an entity in the owned mesh 00547 for( Range::iterator it = owned.begin(); it != owned.end(); it++ ) 00548 { 00549 double volume; 00550 error = vw.quality_measure( *it, q, volume );MB_CHK_ERR( error ); 00551 if( volume > vmax_local ) vmax_local = volume; 00552 } 00553 00554 // Get the global maximum 00555 double vmax_global = vmax_local; 00556 #ifdef MOAB_HAVE_MPI 00557 mpi_err = MPI_Reduce( &vmax_local, &vmax_global, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD ); 00558 if( mpi_err ) 00559 { 00560 MPI_Finalize(); 00561 return MB_FAILURE; 00562 } 00563 #endif 00564 00565 vmax = vmax_global; 00566 00567 return MB_SUCCESS; 00568 } 00569 00570 bool parse_id_list( const char* string, int dim, int nval, std::vector< int >& results ) 00571 { 00572 bool okay = true; 00573 char* mystr = strdup( string ); 00574 for( const char* ptr = strtok( mystr, "," ); ptr; ptr = strtok( 0, "," ) ) 00575 { 00576 char* endptr; 00577 int val = strtol( ptr, &endptr, 0 ); 00578 00579 if( dim == 1 || dim == 2 ) 00580 { 00581 if( val != 2 && val != 3 && val != 5 ) 00582 { 00583 std::cerr << "Not a valid degree for the passed dimension" << std::endl; 00584 okay = false; 00585 break; 00586 } 00587 } 00588 else if( dim == 3 ) 00589 { 00590 if( val != 2 && val != 3 ) 00591 { 00592 std::cerr << "Not a valid degree for the passed dimension" << std::endl; 00593 okay = false; 00594 break; 00595 } 00596 } 00597 00598 if( endptr == ptr || val <= 0 ) 00599 { 00600 std::cerr << "Not a valid id: " << ptr << std::endl; 00601 okay = false; 00602 break; 00603 } 00604 00605 results.push_back( val ); 00606 } 00607 00608 if( (int)results.size() < nval ) 00609 { 00610 for( int i = results.size(); i <= nval - 1; i++ ) 00611 results.push_back( results[0] ); 00612 } 00613 else if( (int)results.size() > nval ) 00614 { 00615 for( int i = results.size(); i > nval; i-- ) 00616 results.pop_back(); 00617 } 00618 00619 free( mystr ); 00620 return okay; 00621 } 00622 00623 bool make_opts_string( std::vector< std::string > options, std::string& opts ) 00624 { 00625 opts.clear(); 00626 if( options.empty() ) return true; 00627 00628 // choose a separator character 00629 std::vector< std::string >::const_iterator i; 00630 char separator = '\0'; 00631 const char* alt_separators = ";+,:\t\n"; 00632 for( const char* sep_ptr = alt_separators; *sep_ptr; ++sep_ptr ) 00633 { 00634 bool seen = false; 00635 for( i = options.begin(); i != options.end(); ++i ) 00636 if( i->find( *sep_ptr, 0 ) != std::string::npos ) 00637 { 00638 seen = true; 00639 break; 00640 } 00641 if( !seen ) 00642 { 00643 separator = *sep_ptr; 00644 break; 00645 } 00646 } 00647 if( !separator ) 00648 { 00649 std::cerr << "Error: cannot find separator character for options string" << std::endl; 00650 return false; 00651 } 00652 if( separator != ';' ) 00653 { 00654 opts = ";"; 00655 opts += separator; 00656 } 00657 00658 // concatenate options 00659 i = options.begin(); 00660 opts += *i; 00661 for( ++i; i != options.end(); ++i ) 00662 { 00663 opts += separator; 00664 opts += *i; 00665 } 00666 00667 return true; 00668 }