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