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