MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }