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