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