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