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