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