MOAB: Mesh Oriented datABase  (version 5.2.1)
umr.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines