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 #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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines