MOAB: Mesh Oriented datABase  (version 5.3.1)
convert.cpp
Go to the documentation of this file.
00001 /**
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 // If Microsoft compiler, then WIN32
00017 #ifndef WIN32
00018 #define WIN32
00019 #endif
00020 
00021 #include "moab/Core.hpp"
00022 #include "moab/Range.hpp"
00023 #include "MBTagConventions.hpp"
00024 #include "moab/ReaderWriterSet.hpp"
00025 #include "moab/ReorderTool.hpp"
00026 #include <iostream>
00027 #include <fstream>
00028 #include <sstream>
00029 #include <iomanip>
00030 #include <set>
00031 #include <list>
00032 #include <cstdlib>
00033 #include <algorithm>
00034 #ifndef WIN32
00035 #include <sys/times.h>
00036 #include <limits.h>
00037 #include <unistd.h>
00038 #endif
00039 #include <ctime>
00040 #ifdef MOAB_HAVE_MPI
00041 #include "moab/ParallelComm.hpp"
00042 #endif
00043 
00044 #ifdef MOAB_HAVE_TEMPESTREMAP
00045 #include "moab/Remapping/TempestRemapper.hpp"
00046 #endif
00047 
00048 #include <cstdio>
00049 
00050 /* Exit values */
00051 #define USAGE_ERROR   1
00052 #define READ_ERROR    2
00053 #define WRITE_ERROR   3
00054 #define OTHER_ERROR   4
00055 #define ENT_NOT_FOUND 5
00056 
00057 using namespace moab;
00058 
00059 static void print_usage( const char* name, std::ostream& stream )
00060 {
00061     stream << "Usage: " << name
00062            << " [-a <sat_file>|-A] [-t] [subset options] [-f format] <input_file> [<input_file2> "
00063               "...] <output_file>"
00064            << std::endl
00065            << "\t-f <format>    - Specify output file format" << std::endl
00066            << "\t-a <acis_file> - ACIS SAT file dumped by .cub reader (same as \"-o "
00067               "SAT_FILE=acis_file\""
00068            << std::endl
00069            << "\t-A             - .cub file reader should not dump a SAT file (depricated default)" << std::endl
00070            << "\t-o option      - Specify write option." << std::endl
00071            << "\t-O option      - Specify read option." << std::endl
00072            << "\t-t             - Time read and write of files." << std::endl
00073            << "\t-g             - Enable verbose/debug output." << std::endl
00074            << "\t-h             - Print this help text and exit." << std::endl
00075            << "\t-l             - List available file formats and exit." << std::endl
00076            << "\t-I <dim>       - Generate internal entities of specified dimension." << std::endl
00077 #ifdef MOAB_HAVE_MPI
00078            << "\t-P             - Append processor ID to output file name" << std::endl
00079            << "\t-p             - Replace '%' with processor ID in input and output file name" << std::endl
00080            << "\t-M[0|1|2]      - Read/write in parallel, optionally also doing "
00081               "resolve_shared_ents (1) and exchange_ghosts (2)"
00082            << std::endl
00083            << "\t-z <file>      - Read metis partition information corresponding to an MPAS grid "
00084               "file and create h5m partition file"
00085            << std::endl
00086 #endif
00087 
00088 #ifdef MOAB_HAVE_TEMPESTREMAP
00089            << "\t-B             - Use TempestRemap exodus file reader and convert to MOAB format" << std::endl
00090            << "\t-b             - Convert MOAB mesh to TempestRemap exodus file writer" << std::endl
00091            << "\t-i             - Name of the global DoF tag to use with mbtempest" << std::endl
00092            << "\t-r             - Order of field DoF (discretization) data; FV=1,SE=[1,N]" << std::endl
00093 #endif
00094            << "\t--             - treat all subsequent options as file names" << std::endl
00095            << "\t                 (allows file names beginning with '-')" << std::endl
00096            << "  subset options: " << std::endl
00097            << "\tEach of the following options should be followed by " << std::endl
00098            << "\ta list of ids.  IDs may be separated with commas.  " << std::endl
00099            << "\tRanges of IDs may be specified with a '-' between " << std::endl
00100            << "\ttwo values.  The list may not contain spaces." << std::endl
00101            << "\t-v  - Volume" << std::endl
00102            << "\t-s  - Surface" << std::endl
00103            << "\t-c  - Curve" << std::endl
00104            << "\t-V  - Vertex" << std::endl
00105            << "\t-m  - Material set (block)" << std::endl
00106            << "\t-d  - Dirichlet set (nodeset)" << std::endl
00107            << "\t-n  - Neumann set (sideset)" << std::endl
00108            << "\t-D  - Parallel partitioning set (PARALLEL_PARTITION)" << std::endl
00109            << "\tThe presence of one or more of the following flags limits " << std::endl
00110            << "\tthe exported mesh to only elements of the corresponding " << std::endl
00111            << "\tdimension.  Vertices are always exported." << std::endl
00112            << "\t-1  - Edges " << std::endl
00113            << "\t-2  - Tri, Quad, Polygon " << std::endl
00114            << "\t-3  - Tet, Hex, Prism, etc. " << std::endl;
00115 }
00116 
00117 static void print_help( const char* name )
00118 {
00119     std::cout << " This program can be used to convert between mesh file\n"
00120                  " formats, extract a subset of a mesh file to a separate\n"
00121                  " file, or both.  The type of file to write is determined\n"
00122                  " from the file extension (e.g. \".vtk\") portion of the\n"
00123                  " output file name.\n"
00124                  " \n"
00125                  " While MOAB strives to export and import all data from\n"
00126                  " each supported file format, most file formats do\n"
00127                  " not support MOAB's entire data model.  Thus MOAB cannot\n"
00128                  " guarantee lossless conversion for any file formats\n"
00129                  " other than the native HDF5 representation.\n"
00130                  "\n";
00131 
00132     print_usage( name, std::cout );
00133     exit( 0 );
00134 }
00135 
00136 static void usage_error( const char* name )
00137 {
00138     print_usage( name, std::cerr );
00139 #ifdef MOAB_HAVE_MPI
00140     MPI_Finalize();
00141 #endif
00142     exit( USAGE_ERROR );
00143 }
00144 
00145 static void list_formats( Interface* );
00146 static bool parse_id_list( const char* string, std::set< int >& );
00147 static void print_id_list( const char*, std::ostream& stream, const std::set< int >& list );
00148 static void reset_times();
00149 static void write_times( std::ostream& stream );
00150 static void remove_entities_from_sets( Interface* gMB, Range& dead_entities, Range& empty_sets );
00151 static void remove_from_vector( std::vector< EntityHandle >& vect, const Range& ents_to_remove );
00152 static bool make_opts_string( std::vector< std::string > options, std::string& result );
00153 static std::string percent_subst( const std::string& s, int val );
00154 
00155 static int process_partition_file( Interface* gMB, std::string& metis_partition_file );
00156 
00157 int main( int argc, char* argv[] )
00158 {
00159     int proc_id = 0;
00160 #ifdef MOAB_HAVE_MPI
00161     MPI_Init( &argc, &argv );
00162     MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
00163 #endif
00164 
00165 #ifdef MOAB_HAVE_TEMPESTREMAP
00166     bool tempestin = false, tempestout = false;
00167 #endif
00168 
00169     Core core;
00170     Interface* gMB = &core;
00171     ErrorCode result;
00172     Range range;
00173 
00174 #if( defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_TEMPESTREMAP ) )
00175     moab::ParallelComm* pcomm = new moab::ParallelComm( gMB, MPI_COMM_WORLD, 0 );
00176 #endif
00177 
00178     bool append_rank        = false;
00179     bool percent_rank_subst = false;
00180     int i, dim;
00181     std::list< std::string >::iterator j;
00182     bool dims[4]       = { false, false, false, false };
00183     const char* format = NULL;    // output file format
00184     std::list< std::string > in;  // input file name list
00185     std::string out;              // output file name
00186     bool verbose = false;
00187     std::set< int > geom[4], mesh[4];      // user-specified IDs
00188     std::vector< EntityHandle > set_list;  // list of user-specified sets to write
00189     std::vector< std::string > write_opts, read_opts;
00190     std::string metis_partition_file;
00191 #ifdef MOAB_HAVE_TEMPESTREMAP
00192     std::string globalid_tag_name;
00193     int spectral_order = 1;
00194 #endif
00195 
00196     const char* const mesh_tag_names[] = { DIRICHLET_SET_TAG_NAME, NEUMANN_SET_TAG_NAME, MATERIAL_SET_TAG_NAME,
00197                                            PARALLEL_PARTITION_TAG_NAME };
00198     const char* const geom_names[]     = { "VERTEX", "CURVE", "SURFACE", "VOLUME" };
00199 
00200     // scan arguments
00201     bool do_flag     = true;
00202     bool print_times = false;
00203     bool generate[]  = { false, false, false };
00204     bool pval;
00205     bool parallel = false, resolve_shared = false, exchange_ghosts = false;
00206     for( i = 1; i < argc; i++ )
00207     {
00208         if( !argv[i][0] ) usage_error( argv[0] );
00209 
00210         if( do_flag && argv[i][0] == '-' )
00211         {
00212             if( !argv[i][1] || ( argv[i][1] != 'M' && argv[i][2] ) ) usage_error( argv[0] );
00213 
00214             switch( argv[i][1] )
00215             {
00216                     // do flag arguments:
00217                 case '-':
00218                     do_flag = false;
00219                     break;
00220                 case 'g':
00221                     verbose = true;
00222                     break;
00223                 case 't':
00224                     print_times = true;
00225                     break;
00226                 case 'A':
00227                     break;
00228                 case 'h':
00229                 case 'H':
00230                     print_help( argv[0] );
00231                     break;
00232                 case 'l':
00233                     list_formats( gMB );
00234                     break;
00235 #ifdef MOAB_HAVE_MPI
00236                 case 'P':
00237                     append_rank = true;
00238                     break;
00239                 case 'p':
00240                     percent_rank_subst = true;
00241                     break;
00242                 case 'M':
00243                     parallel = true;
00244                     if( argv[i][2] == '1' || argv[i][2] == '2' ) resolve_shared = true;
00245                     if( argv[i][2] == '2' ) exchange_ghosts = true;
00246                     break;
00247 #endif
00248 #ifdef MOAB_HAVE_TEMPESTREMAP
00249                 case 'B':
00250                     tempestin = true;
00251                     break;
00252                 case 'b':
00253                     tempestout = true;
00254                     break;
00255 #endif
00256                 case '1':
00257                 case '2':
00258                 case '3':
00259                     dims[argv[i][1] - '0'] = true;
00260                     break;
00261                     // do options that require additional args:
00262                 default:
00263                     ++i;
00264                     if( i == argc || argv[i][0] == '-' )
00265                     {
00266                         std::cerr << "Expected argument following " << argv[i - 1] << std::endl;
00267                         usage_error( argv[0] );
00268                     }
00269                     if( argv[i - 1][1] == 'I' )
00270                     {
00271                         dim = atoi( argv[i] );
00272                         if( dim < 1 || dim > 2 )
00273                         {
00274                             std::cerr << "Invalid dimension value following -I" << std::endl;
00275                             usage_error( argv[0] );
00276                         }
00277                         generate[dim] = true;
00278                         continue;
00279                     }
00280                     pval = false;
00281                     switch( argv[i - 1][1] )
00282                     {
00283                         case 'a':
00284                             read_opts.push_back( std::string( "SAT_FILE=" ) + argv[i] );
00285                             pval = true;
00286                             break;
00287                         case 'f':
00288                             format = argv[i];
00289                             pval   = true;
00290                             break;
00291                         case 'o':
00292                             write_opts.push_back( argv[i] );
00293                             pval = true;
00294                             break;
00295                         case 'O':
00296                             read_opts.push_back( argv[i] );
00297                             pval = true;
00298                             break;
00299 #ifdef MOAB_HAVE_TEMPESTREMAP
00300                         case 'i':
00301                             globalid_tag_name = std::string( argv[i] );
00302                             pval              = true;
00303                             break;
00304                         case 'r':
00305                             spectral_order = atoi( argv[i] );
00306                             pval           = true;
00307                             break;
00308 #endif
00309                         case 'v':
00310                             pval = parse_id_list( argv[i], geom[3] );
00311                             break;
00312                         case 's':
00313                             pval = parse_id_list( argv[i], geom[2] );
00314                             break;
00315                         case 'c':
00316                             pval = parse_id_list( argv[i], geom[1] );
00317                             break;
00318                         case 'V':
00319                             pval = parse_id_list( argv[i], geom[0] );
00320                             break;
00321                         case 'D':
00322                             pval = parse_id_list( argv[i], mesh[3] );
00323                             break;
00324                         case 'm':
00325                             pval = parse_id_list( argv[i], mesh[2] );
00326                             break;
00327                         case 'n':
00328                             pval = parse_id_list( argv[i], mesh[1] );
00329                             break;
00330                         case 'd':
00331                             pval = parse_id_list( argv[i], mesh[0] );
00332                             break;
00333                         case 'z':
00334                             metis_partition_file = argv[i];
00335                             pval                 = true;
00336                             break;
00337                         default:
00338                             std::cerr << "Invalid option: " << argv[i] << std::endl;
00339                     }
00340 
00341                     if( !pval )
00342                     {
00343                         std::cerr << "Invalid flag or flag value: " << argv[i - 1] << " " << argv[i] << std::endl;
00344                         usage_error( argv[0] );
00345                     }
00346             }
00347         }
00348         // do file names
00349         else
00350         {
00351             in.push_back( argv[i] );
00352         }
00353     }
00354     if( in.size() < 2 )
00355     {
00356         std::cerr << "No output file name specified." << std::endl;
00357         usage_error( argv[0] );
00358     }
00359     // output file name is the last one specified
00360     out = in.back();
00361     in.pop_back();
00362 
00363     if( append_rank )
00364     {
00365         std::ostringstream mod;
00366         mod << out << "." << proc_id;
00367         out = mod.str();
00368     }
00369 
00370     if( percent_rank_subst )
00371     {
00372         for( j = in.begin(); j != in.end(); ++j )
00373             *j = percent_subst( *j, proc_id );
00374         out = percent_subst( out, proc_id );
00375     }
00376 
00377     // construct options string from individual options
00378     std::string read_options, write_options;
00379     if( parallel )
00380     {
00381         read_opts.push_back( "PARALLEL=READ_PART" );
00382         read_opts.push_back( "PARTITION=PARALLEL_PARTITION" );
00383         if( !append_rank && !percent_rank_subst ) write_opts.push_back( "PARALLEL=WRITE_PART" );
00384     }
00385     if( resolve_shared ) read_opts.push_back( "PARALLEL_RESOLVE_SHARED_ENTS" );
00386     if( exchange_ghosts ) read_opts.push_back( "PARALLEL_GHOSTS=3.0.1" );
00387 
00388     if( !make_opts_string( read_opts, read_options ) || !make_opts_string( write_opts, write_options ) )
00389     {
00390 #ifdef MOAB_HAVE_MPI
00391         MPI_Finalize();
00392 #endif
00393         return USAGE_ERROR;
00394     }
00395 
00396     if( !metis_partition_file.empty() )
00397     {
00398         if( ( in.size() != 1 ) || ( proc_id != 0 ) )
00399         {
00400             std::cerr << " mpas partition allows only one input file, in serial conversion\n";
00401 #ifdef MOAB_HAVE_MPI
00402             MPI_Finalize();
00403 #endif
00404             return USAGE_ERROR;
00405         }
00406     }
00407 
00408     Tag id_tag = gMB->globalId_tag();
00409 
00410     // Read the input file.
00411 #ifdef MOAB_HAVE_TEMPESTREMAP
00412     if( tempestin && in.size() > 1 )
00413     {
00414         std::cerr << " we can read only one tempest files at a time\n";
00415 #ifdef MOAB_HAVE_MPI
00416         MPI_Finalize();
00417 #endif
00418         return USAGE_ERROR;
00419     }
00420 
00421 #ifdef MOAB_HAVE_MPI
00422     TempestRemapper* remapper = new moab::TempestRemapper( gMB, pcomm );
00423 #else
00424     TempestRemapper* remapper = new moab::TempestRemapper( gMB );
00425 #endif
00426 
00427     bool use_overlap_context = false;
00428     Tag srcParentTag, tgtParentTag;
00429 
00430 #endif
00431     for( j = in.begin(); j != in.end(); ++j )
00432     {
00433         std::string inFileName = *j;
00434 
00435         reset_times();
00436 
00437 #ifdef MOAB_HAVE_TEMPESTREMAP
00438 
00439         remapper->meshValidate = false;
00440         // remapper->constructEdgeMap = true;
00441         remapper->initialize();
00442 
00443         if( tempestin )
00444         {
00445             // convert
00446             result = remapper->LoadMesh( moab::Remapper::SourceMesh, inFileName, moab::TempestRemapper::DEFAULT );MB_CHK_ERR( result );
00447 
00448             Mesh* tempestMesh = remapper->GetMesh( moab::Remapper::SourceMesh );
00449             tempestMesh->RemoveZeroEdges();
00450             tempestMesh->RemoveCoincidentNodes();
00451 
00452             // Load the meshes and validate
00453             result = remapper->ConvertTempestMesh( moab::Remapper::SourceMesh );
00454 
00455             // Check if we are converting a RLL grid
00456             NcFile ncInput( inFileName.c_str(), NcFile::ReadOnly );
00457             bool isRectilinearGrid = false;
00458 
00459             NcError error_temp( NcError::silent_nonfatal );
00460             // get the attribute
00461             NcAtt* attRectilinear = ncInput.get_att( "rectilinear" );
00462 
00463             // If rectilinear attribute present, mark it
00464             std::vector< int > vecDimSizes( 3, 0 );
00465             Tag rectilinearTag;
00466             // Tag data contains: guessed mesh type,     mesh size1,     mesh size 2
00467             //          Example:  CS(0)/ICO(1)/ICOD(2),  num_elements,   num_nodes
00468             //                 :       RLL(3),           num_lat,        num_lon
00469             result = gMB->tag_get_handle( "ClimateMetadata", 3, MB_TYPE_INTEGER, rectilinearTag,
00470                                           MB_TAG_SPARSE | MB_TAG_CREAT, vecDimSizes.data() );MB_CHK_SET_ERR( result, "can't create rectilinear sizes tag" );
00471 
00472             if( attRectilinear != nullptr )
00473             {
00474                 isRectilinearGrid = true;
00475 
00476                 // Obtain rectilinear attributes (dimension sizes)
00477                 NcAtt* attRectilinearDim0Size = ncInput.get_att( "rectilinear_dim0_size" );
00478                 NcAtt* attRectilinearDim1Size = ncInput.get_att( "rectilinear_dim1_size" );
00479 
00480                 if( attRectilinearDim0Size == nullptr )
00481                 {
00482                     _EXCEPTIONT( "Missing attribute \"rectilinear_dim0_size\"" );
00483                 }
00484                 if( attRectilinearDim1Size == nullptr )
00485                 {
00486                     _EXCEPTIONT( "Missing attribute \"rectilinear_dim1_size\"" );
00487                 }
00488 
00489                 int nDim0Size = attRectilinearDim0Size->as_int( 0 );
00490                 int nDim1Size = attRectilinearDim1Size->as_int( 0 );
00491 
00492                 // Obtain rectilinear attributes (dimension names)
00493                 NcAtt* attRectilinearDim0Name = ncInput.get_att( "rectilinear_dim0_name" );
00494                 NcAtt* attRectilinearDim1Name = ncInput.get_att( "rectilinear_dim1_name" );
00495 
00496                 if( attRectilinearDim0Name == nullptr )
00497                 {
00498                     _EXCEPTIONT( "Missing attribute \"rectilinear_dim0_name\"" );
00499                 }
00500                 if( attRectilinearDim1Name == nullptr )
00501                 {
00502                     _EXCEPTIONT( "Missing attribute \"rectilinear_dim1_name\"" );
00503                 }
00504 
00505                 std::string strDim0Name = attRectilinearDim0Name->as_string( 0 );
00506                 std::string strDim1Name = attRectilinearDim1Name->as_string( 0 );
00507 
00508                 std::map< std::string, int > vecDimNameSizes;
00509                 // Push rectilinear attributes into array
00510                 vecDimNameSizes[strDim0Name] = nDim0Size;
00511                 vecDimNameSizes[strDim1Name] = nDim1Size;
00512                 vecDimSizes[0]               = static_cast< int >( moab::TempestRemapper::RLL );
00513                 vecDimSizes[1]               = vecDimNameSizes["lat"];
00514                 vecDimSizes[2]               = vecDimNameSizes["lon"];
00515 
00516                 printf( "Rectilinear RLL mesh size: (lat) %d X (lon) %d\n", vecDimSizes[1], vecDimSizes[2] );
00517 
00518                 moab::EntityHandle mSet = 0;
00519                 // mSet   = remapper->GetMeshSet( moab::Remapper::SourceMesh );
00520                 result = gMB->tag_set_data( rectilinearTag, &mSet, 1, vecDimSizes.data() );MB_CHK_ERR( result );
00521             }
00522             else
00523             {
00524                 const Range& elems = remapper->GetMeshEntities( moab::Remapper::SourceMesh );
00525                 bool isQuads       = elems.all_of_type( moab::MBQUAD );
00526                 bool isTris        = elems.all_of_type( moab::MBTRI );
00527                 // vecDimSizes[0] = static_cast< int >( remapper->GetMeshType( moab::Remapper::SourceMesh );
00528                 vecDimSizes[0] = ( isQuads ? static_cast< int >( moab::TempestRemapper::CS )
00529                                            : ( isTris ? static_cast< int >( moab::TempestRemapper::ICO )
00530                                                       : static_cast< int >( moab::TempestRemapper::ICOD ) ) );
00531                 vecDimSizes[1] = elems.size();
00532                 vecDimSizes[2] = remapper->GetMeshVertices( moab::Remapper::SourceMesh ).size();
00533 
00534                 switch( vecDimSizes[0] )
00535                 {
00536                     case 0:
00537                         printf( "Cubed-Sphere mesh: %d (elems), %d (nodes)\n", vecDimSizes[1], vecDimSizes[2] );
00538                         break;
00539                     case 2:
00540                         printf( "Icosahedral (triangular) mesh: %d (elems), %d (nodes)\n", vecDimSizes[1],
00541                                 vecDimSizes[2] );
00542                         break;
00543                     case 3:
00544                     default:
00545                         printf( "Polygonal mesh: %d (elems), %d (nodes)\n", vecDimSizes[1], vecDimSizes[2] );
00546                         break;
00547                 }
00548 
00549                 moab::EntityHandle mSet = 0;
00550                 // mSet   = remapper->GetMeshSet( moab::Remapper::SourceMesh );
00551                 result = gMB->tag_set_data( rectilinearTag, &mSet, 1, vecDimSizes.data() );MB_CHK_ERR( result );
00552             }
00553 
00554             ncInput.close();
00555 
00556             const size_t nOverlapFaces = tempestMesh->faces.size();
00557             if( tempestMesh->vecSourceFaceIx.size() == nOverlapFaces &&
00558                 tempestMesh->vecSourceFaceIx.size() == nOverlapFaces )
00559             {
00560                 int defaultInt      = -1;
00561                 use_overlap_context = true;
00562                 // Check if our MOAB mesh has RED and BLUE tags; this would indicate we are
00563                 // converting an overlap grid
00564                 result = gMB->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag,
00565                                               MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( result, "can't create target parent tag" );
00566 
00567                 result = gMB->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag,
00568                                               MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( result, "can't create source parent tag" );
00569 
00570                 const Range& faces = remapper->GetMeshEntities( moab::Remapper::SourceMesh );
00571 
00572                 std::vector< int > gids( faces.size() ), srcpar( faces.size() ), tgtpar( faces.size() );
00573                 result = gMB->tag_get_data( id_tag, faces, &gids[0] );MB_CHK_ERR( result );
00574 
00575                 for( unsigned ii = 0; ii < faces.size(); ++ii )
00576                 {
00577                     srcpar[ii] = tempestMesh->vecSourceFaceIx[gids[ii] - 1];
00578                     tgtpar[ii] = tempestMesh->vecTargetFaceIx[gids[ii] - 1];
00579                 }
00580 
00581                 result = gMB->tag_set_data( srcParentTag, faces, &srcpar[0] );MB_CHK_ERR( result );
00582                 result = gMB->tag_set_data( tgtParentTag, faces, &tgtpar[0] );MB_CHK_ERR( result );
00583 
00584                 srcpar.clear();
00585                 tgtpar.clear();
00586                 gids.clear();
00587             }
00588         }
00589         else if( tempestout )
00590         {
00591             moab::EntityHandle& srcmesh = remapper->GetMeshSet( moab::Remapper::SourceMesh );
00592             moab::EntityHandle& ovmesh  = remapper->GetMeshSet( moab::Remapper::OverlapMesh );
00593 
00594             // load the mesh in MOAB format
00595             std::vector< int > metadata;
00596             result = remapper->LoadNativeMesh( *j, srcmesh, metadata );MB_CHK_ERR( result );
00597 
00598             // Check if our MOAB mesh has RED and BLUE tags; this would indicate we are converting
00599             // an overlap grid
00600             ErrorCode rval1 = gMB->tag_get_handle( "SourceParent", srcParentTag );
00601             ErrorCode rval2 = gMB->tag_get_handle( "TargetParent", tgtParentTag );
00602             if( rval1 == MB_SUCCESS && rval2 == MB_SUCCESS )
00603             {
00604                 use_overlap_context = true;
00605                 ovmesh              = srcmesh;
00606 
00607                 Tag countTag;
00608                 result = gMB->tag_get_handle( "Counting", countTag );
00609                 // std::vector<int> count_ids()
00610 
00611                 // Load the meshes and validate
00612                 Tag order;
00613                 ReorderTool reorder_tool( &core );
00614                 result = reorder_tool.handle_order_from_int_tag( srcParentTag, -1, order );MB_CHK_ERR( result );
00615                 result = reorder_tool.reorder_entities( order );MB_CHK_ERR( result );
00616                 result = gMB->tag_delete( order );MB_CHK_ERR( result );
00617                 result = remapper->ConvertMeshToTempest( moab::Remapper::OverlapMesh );MB_CHK_ERR( result );
00618             }
00619             else
00620             {
00621                 if( metadata[0] == static_cast< int >( moab::TempestRemapper::RLL ) )
00622                 {
00623                     assert( metadata.size() );
00624                     std::cout << "Converting a RLL mesh with rectilinear dimension: " << metadata[0] << " X "
00625                               << metadata[1] << std::endl;
00626                 }
00627 
00628                 // Convert the mesh and validate
00629                 result = remapper->ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( result );
00630             }
00631         }
00632         else
00633             result = gMB->load_file( j->c_str(), 0, read_options.c_str() );
00634 #else
00635         result = gMB->load_file( j->c_str(), 0, read_options.c_str() );
00636 #endif
00637         if( MB_SUCCESS != result )
00638         {
00639             std::cerr << "Failed to load \"" << *j << "\"." << std::endl;
00640             std::cerr << "Error code: " << gMB->get_error_string( result ) << " (" << result << ")" << std::endl;
00641             std::string message;
00642             if( MB_SUCCESS == gMB->get_last_error( message ) && !message.empty() )
00643                 std::cerr << "Error message: " << message << std::endl;
00644 #ifdef MOAB_HAVE_MPI
00645             MPI_Finalize();
00646 #endif
00647             return READ_ERROR;
00648         }
00649         if( !proc_id ) std::cerr << "Read \"" << *j << "\"" << std::endl;
00650         if( print_times && !proc_id ) write_times( std::cout );
00651     }
00652 
00653     // Determine if the user has specified any geometry sets to write
00654     bool have_geom = false;
00655     for( dim = 0; dim <= 3; ++dim )
00656     {
00657         if( !geom[dim].empty() ) have_geom = true;
00658         if( verbose ) print_id_list( geom_names[dim], std::cout, geom[dim] );
00659     }
00660 
00661     // True if the user has specified any sets to write
00662     bool have_sets = have_geom;
00663 
00664     // Get geometry tags
00665     Tag dim_tag;
00666     if( have_geom )
00667     {
00668         if( id_tag == 0 )
00669         {
00670             std::cerr << "No ID tag defined." << std::endl;
00671             have_geom = false;
00672         }
00673         result = gMB->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dim_tag );
00674         if( MB_SUCCESS != result )
00675         {
00676             std::cerr << "No geometry tag defined." << std::endl;
00677             have_geom = false;
00678         }
00679     }
00680 
00681     // Get geometry sets
00682     if( have_geom )
00683     {
00684         int id_val;
00685         Tag tags[]         = { id_tag, dim_tag };
00686         const void* vals[] = { &id_val, &dim };
00687         for( dim = 0; dim <= 3; ++dim )
00688         {
00689             int init_count = set_list.size();
00690             for( std::set< int >::iterator iter = geom[dim].begin(); iter != geom[dim].end(); ++iter )
00691             {
00692                 id_val = *iter;
00693                 range.clear();
00694                 result = gMB->get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, range );
00695                 if( MB_SUCCESS != result || range.empty() )
00696                 {
00697                     range.clear();
00698                     std::cerr << geom_names[dim] << " " << id_val << " not found.\n";
00699                 }
00700                 std::copy( range.begin(), range.end(), std::back_inserter( set_list ) );
00701             }
00702 
00703             if( verbose )
00704                 std::cout << "Found " << ( set_list.size() - init_count ) << ' ' << geom_names[dim] << " sets"
00705                           << std::endl;
00706         }
00707     }
00708 
00709     // Get mesh groupings
00710     for( i = 0; i < 4; ++i )
00711     {
00712         if( verbose ) print_id_list( mesh_tag_names[i], std::cout, mesh[i] );
00713 
00714         if( mesh[i].empty() ) continue;
00715         have_sets = true;
00716 
00717         // Get tag
00718         Tag tag;
00719         result = gMB->tag_get_handle( mesh_tag_names[i], 1, MB_TYPE_INTEGER, tag );
00720         if( MB_SUCCESS != result )
00721         {
00722             std::cerr << "Tag not found: " << mesh_tag_names[i] << std::endl;
00723             continue;
00724         }
00725 
00726         // get entity sets
00727         int init_count = set_list.size();
00728         for( std::set< int >::iterator iter = mesh[i].begin(); iter != mesh[i].end(); ++iter )
00729         {
00730             range.clear();
00731             const void* vals[] = { &*iter };
00732             result             = gMB->get_entities_by_type_and_tag( 0, MBENTITYSET, &tag, vals, 1, range );
00733             if( MB_SUCCESS != result || range.empty() )
00734             {
00735                 range.clear();
00736                 std::cerr << mesh_tag_names[i] << " " << *iter << " not found.\n";
00737             }
00738             std::copy( range.begin(), range.end(), std::back_inserter( set_list ) );
00739         }
00740 
00741         if( verbose )
00742             std::cout << "Found " << ( set_list.size() - init_count ) << ' ' << mesh_tag_names[i] << " sets"
00743                       << std::endl;
00744     }
00745 
00746     // Check if output is limited to certain dimensions of elements
00747     bool bydim = false;
00748     for( dim = 1; dim < 4; ++dim )
00749         if( dims[dim] ) bydim = true;
00750 
00751     // Check conflicting input
00752     if( bydim )
00753     {
00754         if( generate[1] && !dims[1] )
00755         {
00756             std::cerr << "Warning: Request to generate 1D internal entities but not export them." << std::endl;
00757             generate[1] = false;
00758         }
00759         if( generate[2] && !dims[2] )
00760         {
00761             std::cerr << "Warning: Request to generate 2D internal entities but not export them." << std::endl;
00762             generate[2] = false;
00763         }
00764     }
00765 
00766     // Generate any internal entities
00767     if( generate[1] || generate[2] )
00768     {
00769         EntityHandle all_mesh    = 0;
00770         const EntityHandle* sets = &all_mesh;
00771         int num_sets             = 1;
00772         if( have_sets )
00773         {
00774             num_sets = set_list.size();
00775             sets     = &set_list[0];
00776         }
00777         for( i = 0; i < num_sets; ++i )
00778         {
00779             Range dim3, dim2, adj;
00780             gMB->get_entities_by_dimension( sets[i], 3, dim3, true );
00781             if( generate[1] )
00782             {
00783                 gMB->get_entities_by_dimension( sets[i], 2, dim2, true );
00784                 gMB->get_adjacencies( dim3, 1, true, adj, Interface::UNION );
00785                 gMB->get_adjacencies( dim2, 1, true, adj, Interface::UNION );
00786             }
00787             if( generate[2] )
00788             {
00789                 gMB->get_adjacencies( dim3, 2, true, adj, Interface::UNION );
00790             }
00791             if( sets[i] ) gMB->add_entities( sets[i], adj );
00792         }
00793     }
00794 
00795     // Delete any entities not of the dimensions to be exported
00796     if( bydim )
00797     {
00798         // Get list of dead elements
00799         Range dead_entities, tmp_range;
00800         for( dim = 1; dim <= 3; ++dim )
00801         {
00802             if( dims[dim] ) continue;
00803             gMB->get_entities_by_dimension( 0, dim, tmp_range );
00804             dead_entities.merge( tmp_range );
00805         }
00806         // Remove dead entities from all sets, and add all
00807         // empty sets to list of dead entities.
00808         Range empty_sets;
00809         remove_entities_from_sets( gMB, dead_entities, empty_sets );
00810         while( !empty_sets.empty() )
00811         {
00812             if( !set_list.empty() ) remove_from_vector( set_list, empty_sets );
00813             dead_entities.merge( empty_sets );
00814             tmp_range.clear();
00815             remove_entities_from_sets( gMB, empty_sets, tmp_range );
00816             empty_sets = subtract( tmp_range, dead_entities );
00817         }
00818         // Destroy dead entities
00819         gMB->delete_entities( dead_entities );
00820     }
00821 
00822     // If user specified sets to write, but none were found, exit.
00823     if( have_sets && set_list.empty() )
00824     {
00825         std::cerr << "Nothing to write." << std::endl;
00826 #ifdef MOAB_HAVE_MPI
00827         MPI_Finalize();
00828 #endif
00829         return ENT_NOT_FOUND;
00830     }
00831 
00832     // interpret the mpas partition file created by gpmetis
00833     if( !metis_partition_file.empty() )
00834     {
00835         int err = process_partition_file( gMB, metis_partition_file );
00836         if( err )
00837         {
00838             std::cerr << "Failed to process partition file \"" << metis_partition_file << "\"." << std::endl;
00839 #ifdef MOAB_HAVE_MPI
00840             MPI_Finalize();
00841 #endif
00842             return WRITE_ERROR;
00843         }
00844     }
00845     if( verbose )
00846     {
00847         if( have_sets )
00848             std::cout << "Found " << set_list.size() << " specified sets to write (total)." << std::endl;
00849         else
00850             std::cout << "No sets specifed.  Writing entire mesh." << std::endl;
00851     }
00852 
00853     // Write the output file
00854     reset_times();
00855 #ifdef MOAB_HAVE_TEMPESTREMAP
00856     Range faces;
00857     Mesh* tempestMesh =
00858         remapper->GetMesh( ( use_overlap_context ? moab::Remapper::OverlapMesh : moab::Remapper::SourceMesh ) );
00859     moab::EntityHandle& srcmesh =
00860         remapper->GetMeshSet( ( use_overlap_context ? moab::Remapper::OverlapMesh : moab::Remapper::SourceMesh ) );
00861     result = gMB->get_entities_by_dimension( srcmesh, 2, faces );MB_CHK_ERR( result );
00862     int ntot_elements = 0, nelements = faces.size();
00863 #ifdef MOAB_HAVE_MPI
00864     int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, pcomm->comm() );
00865     if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
00866 #else
00867     ntot_elements             = nelements;
00868 #endif
00869 
00870     Tag gidTag = gMB->globalId_tag();
00871     std::vector< int > gids( faces.size() );
00872     result = gMB->tag_get_data( gidTag, faces, &gids[0] );MB_CHK_ERR( result );
00873 
00874     if( faces.size() > 1 && gids[0] == gids[1] && !use_overlap_context )
00875     {
00876 #ifdef MOAB_HAVE_MPI
00877         result = pcomm->assign_global_ids( srcmesh, 2, 1, false );MB_CHK_ERR( result );
00878 #else
00879         result = remapper->assign_vertex_element_IDs( gidTag, srcmesh, 2, 1 );MB_CHK_ERR( result );
00880         result = remapper->assign_vertex_element_IDs( gidTag, srcmesh, 0, 1 );MB_CHK_ERR( result );
00881 #endif
00882     }
00883 
00884     // VSM: If user requested explicitly for some metadata, we need to generate the DoF ID tag
00885     // and set the appropriate numbering based on specified discretization order
00886     // Useful only for SE meshes with GLL DoFs
00887     if( spectral_order > 1 && globalid_tag_name.size() > 1 )
00888     {
00889         result = remapper->GenerateMeshMetadata( *tempestMesh, ntot_elements, faces, NULL, globalid_tag_name,
00890                                                  spectral_order );MB_CHK_ERR( result );
00891     }
00892 
00893     if( tempestout )
00894     {
00895         // Check if our MOAB mesh has RED and BLUE tags; this would indicate we are converting an
00896         // overlap grid
00897         if( use_overlap_context && false )
00898         {
00899             const int nOverlapFaces = faces.size();
00900             // Overlap mesh: resize the source and target connection arrays
00901             tempestMesh->vecSourceFaceIx.resize( nOverlapFaces );  // 0-based indices corresponding to source mesh
00902             tempestMesh->vecTargetFaceIx.resize( nOverlapFaces );  // 0-based indices corresponding to target mesh
00903             result = gMB->tag_get_data( srcParentTag, faces, &tempestMesh->vecSourceFaceIx[0] );MB_CHK_ERR( result );
00904             result = gMB->tag_get_data( tgtParentTag, faces, &tempestMesh->vecTargetFaceIx[0] );MB_CHK_ERR( result );
00905         }
00906         // Write out the mesh using TempestRemap
00907         tempestMesh->Write( out, NcFile::Netcdf4 );
00908     }
00909     else
00910     {
00911 #endif
00912 
00913         if( have_sets )
00914             result = gMB->write_file( out.c_str(), format, write_options.c_str(), &set_list[0], set_list.size() );
00915         else
00916             result = gMB->write_file( out.c_str(), format, write_options.c_str() );
00917         if( MB_SUCCESS != result )
00918         {
00919             std::cerr << "Failed to write \"" << out << "\"." << std::endl;
00920             std::cerr << "Error code: " << gMB->get_error_string( result ) << " (" << result << ")" << std::endl;
00921             std::string message;
00922             if( MB_SUCCESS == gMB->get_last_error( message ) && !message.empty() )
00923                 std::cerr << "Error message: " << message << std::endl;
00924 #ifdef MOAB_HAVE_MPI
00925             MPI_Finalize();
00926 #endif
00927             return WRITE_ERROR;
00928         }
00929 #ifdef MOAB_HAVE_TEMPESTREMAP
00930     }
00931     delete remapper;
00932 #endif
00933 
00934     if( !proc_id ) std::cerr << "Wrote \"" << out << "\"" << std::endl;
00935     if( print_times && !proc_id ) write_times( std::cout );
00936 
00937 #ifdef MOAB_HAVE_MPI
00938     MPI_Finalize();
00939 #endif
00940     return 0;
00941 }
00942 
00943 bool parse_id_list( const char* string, std::set< int >& results )
00944 {
00945     bool okay   = true;
00946     char* mystr = strdup( string );
00947     for( const char* ptr = strtok( mystr, "," ); ptr; ptr = strtok( 0, "," ) )
00948     {
00949         char* endptr;
00950         long val = strtol( ptr, &endptr, 0 );
00951         if( endptr == ptr || val <= 0 )
00952         {
00953             std::cerr << "Not a valid id: " << ptr << std::endl;
00954             okay = false;
00955             break;
00956         }
00957 
00958         long val2 = val;
00959         if( *endptr == '-' )
00960         {
00961             const char* sptr = endptr + 1;
00962             val2             = strtol( sptr, &endptr, 0 );
00963             if( endptr == sptr || val2 <= 0 )
00964             {
00965                 std::cerr << "Not a valid id: " << sptr << std::endl;
00966                 okay = false;
00967                 break;
00968             }
00969             if( val2 < val )
00970             {
00971                 std::cerr << "Invalid id range: " << ptr << std::endl;
00972                 okay = false;
00973                 break;
00974             }
00975         }
00976 
00977         if( *endptr )
00978         {
00979             std::cerr << "Unexpected character: " << *endptr << std::endl;
00980             okay = false;
00981             break;
00982         }
00983 
00984         for( ; val <= val2; ++val )
00985             if( !results.insert( (int)val ).second ) std::cerr << "Warning: duplicate Id: " << val << std::endl;
00986     }
00987 
00988     free( mystr );
00989     return okay;
00990 }
00991 
00992 void print_id_list( const char* head, std::ostream& stream, const std::set< int >& list )
00993 {
00994     stream << head << ": ";
00995 
00996     if( list.empty() )
00997     {
00998         stream << "(none)" << std::endl;
00999         return;
01000     }
01001 
01002     int start, prev;
01003     std::set< int >::const_iterator iter = list.begin();
01004     start = prev = *( iter++ );
01005     for( ;; )
01006     {
01007         if( iter == list.end() || *iter != 1 + prev )
01008         {
01009             stream << start;
01010             if( prev != start ) stream << '-' << prev;
01011             if( iter == list.end() ) break;
01012             stream << ", ";
01013             start = *iter;
01014         }
01015         prev = *( iter++ );
01016     }
01017 
01018     stream << std::endl;
01019 }
01020 
01021 static void print_time( int clk_per_sec, const char* prefix, clock_t ticks, std::ostream& stream )
01022 {
01023     ticks *= clk_per_sec / 100;
01024     clock_t centi   = ticks % 100;
01025     clock_t seconds = ticks / 100;
01026     stream << prefix;
01027     if( seconds < 120 )
01028     {
01029         stream << ( ticks / 100 ) << "." << centi << "s" << std::endl;
01030     }
01031     else
01032     {
01033         clock_t minutes = ( seconds / 60 ) % 60;
01034         clock_t hours   = ( seconds / 3600 );
01035         seconds %= 60;
01036         if( hours ) stream << hours << "h";
01037         if( minutes ) stream << minutes << "m";
01038         if( seconds || centi ) stream << seconds << "." << centi << "s";
01039         stream << " (" << ( ticks / 100 ) << "." << centi << "s)" << std::endl;
01040     }
01041 }
01042 
01043 clock_t usr_time, sys_time, abs_time;
01044 
01045 #ifdef WIN32
01046 
01047 void reset_times()
01048 {
01049     abs_time = clock();
01050 }
01051 
01052 void write_times( std::ostream& stream )
01053 {
01054     clock_t abs_tm = clock();
01055     print_time( CLOCKS_PER_SEC, "  ", abs_tm - abs_time, stream );
01056     abs_time = abs_tm;
01057 }
01058 
01059 #else
01060 
01061 void reset_times()
01062 {
01063     tms timebuf;
01064     abs_time = times( &timebuf );
01065     usr_time = timebuf.tms_utime;
01066     sys_time = timebuf.tms_stime;
01067 }
01068 
01069 void write_times( std::ostream& stream )
01070 {
01071     clock_t usr_tm, sys_tm, abs_tm;
01072     tms timebuf;
01073     abs_tm = times( &timebuf );
01074     usr_tm = timebuf.tms_utime;
01075     sys_tm = timebuf.tms_stime;
01076     print_time( sysconf( _SC_CLK_TCK ), "  real:   ", abs_tm - abs_time, stream );
01077     print_time( sysconf( _SC_CLK_TCK ), "  user:   ", usr_tm - usr_time, stream );
01078     print_time( sysconf( _SC_CLK_TCK ), "  system: ", sys_tm - sys_time, stream );
01079     abs_time = abs_tm;
01080     usr_time = usr_tm;
01081     sys_time = sys_tm;
01082 }
01083 
01084 #endif
01085 
01086 bool make_opts_string( std::vector< std::string > options, std::string& opts )
01087 {
01088     opts.clear();
01089     if( options.empty() ) return true;
01090 
01091     // choose a separator character
01092     std::vector< std::string >::const_iterator i;
01093     char separator             = '\0';
01094     const char* alt_separators = ";+,:\t\n";
01095     for( const char* sep_ptr = alt_separators; *sep_ptr; ++sep_ptr )
01096     {
01097         bool seen = false;
01098         for( i = options.begin(); i != options.end(); ++i )
01099             if( i->find( *sep_ptr, 0 ) != std::string::npos )
01100             {
01101                 seen = true;
01102                 break;
01103             }
01104         if( !seen )
01105         {
01106             separator = *sep_ptr;
01107             break;
01108         }
01109     }
01110     if( !separator )
01111     {
01112         std::cerr << "Error: cannot find separator character for options string" << std::endl;
01113         return false;
01114     }
01115     if( separator != ';' )
01116     {
01117         opts = ";";
01118         opts += separator;
01119     }
01120 
01121     // concatenate options
01122     i = options.begin();
01123     opts += *i;
01124     for( ++i; i != options.end(); ++i )
01125     {
01126         opts += separator;
01127         opts += *i;
01128     }
01129 
01130     return true;
01131 }
01132 
01133 void list_formats( Interface* gMB )
01134 {
01135     const char iface_name[] = "ReaderWriterSet";
01136     ErrorCode err;
01137     ReaderWriterSet* set = 0;
01138     ReaderWriterSet::iterator i;
01139     std::ostream& str = std::cout;
01140 
01141     // get ReaderWriterSet
01142     err = gMB->query_interface( set );
01143     if( err != MB_SUCCESS || !set )
01144     {
01145         std::cerr << "Internal error:  Interface \"" << iface_name << "\" not available.\n";
01146         exit( OTHER_ERROR );
01147     }
01148 
01149     // get field width for format description
01150     size_t w = 0;
01151     for( i = set->begin(); i != set->end(); ++i )
01152         if( i->description().length() > w ) w = i->description().length();
01153 
01154     // write table header
01155     str << "Format  " << std::setw( w ) << std::left << "Description"
01156         << "  Read  Write  File Name Suffixes\n"
01157         << "------  " << std::setw( w ) << std::setfill( '-' ) << "" << std::setfill( ' ' )
01158         << "  ----  -----  ------------------\n";
01159 
01160     // write table data
01161     for( i = set->begin(); i != set->end(); ++i )
01162     {
01163         std::vector< std::string > ext;
01164         i->get_extensions( ext );
01165         str << std::setw( 6 ) << i->name() << "  " << std::setw( w ) << std::left << i->description() << "  "
01166             << ( i->have_reader() ? " yes" : "  no" ) << "  " << ( i->have_writer() ? "  yes" : "   no" ) << " ";
01167         for( std::vector< std::string >::iterator j = ext.begin(); j != ext.end(); ++j )
01168             str << " " << *j;
01169         str << std::endl;
01170     }
01171     str << std::endl;
01172 
01173     gMB->release_interface( set );
01174     exit( 0 );
01175 }
01176 
01177 void remove_entities_from_sets( Interface* gMB, Range& dead_entities, Range& empty_sets )
01178 {
01179     empty_sets.clear();
01180     Range sets;
01181     gMB->get_entities_by_type( 0, MBENTITYSET, sets );
01182     for( Range::iterator i = sets.begin(); i != sets.end(); ++i )
01183     {
01184         Range set_contents;
01185         gMB->get_entities_by_handle( *i, set_contents, false );
01186         set_contents = intersect( set_contents, dead_entities );
01187         gMB->remove_entities( *i, set_contents );
01188         set_contents.clear();
01189         gMB->get_entities_by_handle( *i, set_contents, false );
01190         if( set_contents.empty() ) empty_sets.insert( *i );
01191     }
01192 }
01193 
01194 void remove_from_vector( std::vector< EntityHandle >& vect, const Range& ents_to_remove )
01195 {
01196     Range::const_iterator i;
01197     std::vector< EntityHandle >::iterator j;
01198     for( i = ents_to_remove.begin(); i != ents_to_remove.end(); ++i )
01199     {
01200         j = std::find( vect.begin(), vect.end(), *i );
01201         if( j != vect.end() ) vect.erase( j );
01202     }
01203 }
01204 
01205 std::string percent_subst( const std::string& s, int val )
01206 {
01207     if( s.empty() ) return s;
01208 
01209     size_t j = s.find( '%' );
01210     if( j == std::string::npos ) return s;
01211 
01212     std::ostringstream st;
01213     st << s.substr( 0, j );
01214     st << val;
01215 
01216     size_t i;
01217     while( ( i = s.find( '%', j + 1 ) ) != std::string::npos )
01218     {
01219         st << s.substr( j, i - j );
01220         st << val;
01221         j = i;
01222     }
01223     st << s.substr( j + 1 );
01224     return st.str();
01225 }
01226 
01227 int process_partition_file( Interface* mb, std::string& metis_partition_file )
01228 {
01229     // how many faces in the file ? how do we make sure it is an mpas file?
01230     // mpas atmosphere files can be downloaded from here
01231     // https://mpas-dev.github.io/atmosphere/atmosphere_meshes.html
01232     Range faces;
01233     ErrorCode rval = mb->get_entities_by_dimension( 0, 2, faces );MB_CHK_ERR( rval );
01234     std::cout << " MPAS model has " << faces.size() << " polygons\n";
01235 
01236     // read the partition file
01237     std::ifstream partfile;
01238     partfile.open( metis_partition_file.c_str() );
01239     std::string line;
01240     std::vector< int > parts;
01241     parts.resize( faces.size(), -1 );
01242     int i = 0;
01243     if( partfile.is_open() )
01244     {
01245         while( getline( partfile, line ) )
01246         {
01247             // cout << line << '\n';
01248             parts[i++] = atoi( line.c_str() );
01249             if( i > (int)faces.size() )
01250             {
01251                 std::cerr << " too many lines in partition file \n. bail out \n";
01252                 return 1;
01253             }
01254         }
01255         partfile.close();
01256     }
01257     std::vector< int >::iterator pmax = max_element( parts.begin(), parts.end() );
01258     std::vector< int >::iterator pmin = min_element( parts.begin(), parts.end() );
01259     if( *pmin <= -1 )
01260     {
01261         std::cerr << " partition file is incomplete, *pmin is -1 !! \n";
01262         return 1;
01263     }
01264     std::cout << " partitions range: " << *pmin << " " << *pmax << "\n";
01265     Tag part_set_tag;
01266     int dum_id = -1;
01267     rval = mb->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_set_tag, MB_TAG_SPARSE | MB_TAG_CREAT,
01268                                &dum_id );MB_CHK_ERR( rval );
01269 
01270     // get any sets already with this tag, and clear them
01271     // remove the parallel partition sets if they exist
01272     Range tagged_sets;
01273     rval = mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &part_set_tag, NULL, 1, tagged_sets, Interface::UNION );MB_CHK_ERR( rval );
01274     if( !tagged_sets.empty() )
01275     {
01276         rval = mb->clear_meshset( tagged_sets );MB_CHK_ERR( rval );
01277         rval = mb->tag_delete_data( part_set_tag, tagged_sets );MB_CHK_ERR( rval );
01278     }
01279     Tag gid;
01280     rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_ERR( rval );
01281     int num_sets = *pmax + 1;
01282     if( *pmin != 0 )
01283     {
01284         std::cout << " problem reading parts; min is not 0 \n";
01285         return 1;
01286     }
01287     for( i = 0; i < num_sets; i++ )
01288     {
01289         EntityHandle new_set;
01290         rval = mb->create_meshset( MESHSET_SET, new_set );MB_CHK_ERR( rval );
01291         tagged_sets.insert( new_set );
01292     }
01293     int* dum_ids = new int[num_sets];
01294     for( i = 0; i < num_sets; i++ )
01295         dum_ids[i] = i;
01296 
01297     rval = mb->tag_set_data( part_set_tag, tagged_sets, dum_ids );MB_CHK_ERR( rval );
01298     delete[] dum_ids;
01299 
01300     std::vector< int > gids;
01301     int num_faces = (int)faces.size();
01302     gids.resize( num_faces );
01303     rval = mb->tag_get_data( gid, faces, &gids[0] );MB_CHK_ERR( rval );
01304 
01305     for( int j = 0; j < num_faces; j++ )
01306     {
01307         int eid         = gids[j];
01308         EntityHandle eh = faces[j];
01309         int partition   = parts[eid - 1];
01310         if( partition < 0 || partition >= num_sets )
01311         {
01312             std::cout << " wrong partition number \n";
01313             return 1;
01314         }
01315         rval = mb->add_entities( tagged_sets[partition], &eh, 1 );MB_CHK_ERR( rval );
01316     }
01317     return 0;
01318 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines