 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /**
2 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain rights in this software.
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 */
15
16 // If Microsoft compiler, then WIN32
17 #ifndef WIN32
18 #define WIN32
19 #endif
20
21 #include "moab/Core.hpp"
22 #include "moab/Range.hpp"
23 #include "MBTagConventions.hpp"
24 #include "moab/ReaderWriterSet.hpp"
25 #include "moab/ReorderTool.hpp"
26 #include <iostream>
27 #include <fstream>
28 #include <sstream>
29 #include <iomanip>
30 #include <set>
31 #include <list>
32 #include <cstdlib>
33 #include <algorithm>
34 #ifndef WIN32
35 #include <sys/times.h>
36 #include <limits.h>
37 #include <unistd.h>
38 #endif
39 #include <ctime>
40 #ifdef MOAB_HAVE_MPI
41 #include "moab/ParallelComm.hpp"
42 #endif
43
44 #ifdef MOAB_HAVE_TEMPESTREMAP
45 #include "moab/Remapping/TempestRemapper.hpp"
46
47 constexpr bool offlineGenerator = true;
48 #endif
49
50 #include <cstdio>
51
52 /* Exit values */
53 #define USAGE_ERROR 1
54 #define READ_ERROR 2
55 #define WRITE_ERROR 3
56 #define OTHER_ERROR 4
57 #define ENT_NOT_FOUND 5
58
59 using namespace moab;
60
61 static void print_usage( const char* name, std::ostream& stream )
62 {
63 stream << "Usage: " << name
64 << " [-a <sat_file>|-A] [-t] [subset options] [-f format] <input_file> [<input_file2> "
65 "...] <output_file>"
66 << std::endl
67 << "\t-f <format> - Specify output file format" << std::endl
68 << "\t-a <acis_file> - ACIS SAT file dumped by .cub reader (same as \"-o "
69 "SAT_FILE=acis_file\""
70 << std::endl
71 << "\t-A - .cub file reader should not dump a SAT file (depricated default)" << std::endl
72 << "\t-o option - Specify write option." << std::endl
73 << "\t-O option - Specify read option." << std::endl
74 << "\t-t - Time read and write of files." << std::endl
75 << "\t-g - Enable verbose/debug output." << std::endl
76 << "\t-h - Print this help text and exit." << std::endl
77 << "\t-l - List available file formats and exit." << std::endl
78 << "\t-I <dim> - Generate internal entities of specified dimension." << std::endl
79 #ifdef MOAB_HAVE_MPI
80 << "\t-P - Append processor ID to output file name" << std::endl
81 << "\t-p - Replace '%' with processor ID in input and output file name" << std::endl
82 << "\t-M[0|1|2] - Read/write in parallel, optionally also doing "
83 "resolve_shared_ents (1) and exchange_ghosts (2)"
84 << std::endl
85 << "\t-z <file> - Read metis partition information corresponding to an MPAS grid "
86 "file and create h5m partition file"
87 << std::endl
88 #endif
89
90 #ifdef MOAB_HAVE_TEMPESTREMAP
91 << "\t-B - Use TempestRemap exodus file reader and convert to MOAB format" << std::endl
92 << "\t-b - Convert MOAB mesh to TempestRemap exodus file writer" << std::endl
93 << "\t-S - Scale climate mesh to unit sphere" << std::endl
94 << "\t-i - Name of the global DoF tag to use with mbtempest" << std::endl
95 << "\t-r - Order of field DoF (discretization) data; FV=1,SE=[1,N]" << std::endl
96 #endif
97 << "\t-- - treat all subsequent options as file names" << std::endl
98 << "\t (allows file names beginning with '-')" << std::endl
99 << " subset options: " << std::endl
100 << "\tEach of the following options should be followed by " << std::endl
101 << "\ta list of ids. IDs may be separated with commas. " << std::endl
102 << "\tRanges of IDs may be specified with a '-' between " << std::endl
103 << "\ttwo values. The list may not contain spaces." << std::endl
104 << "\t-v - Volume" << std::endl
105 << "\t-s - Surface" << std::endl
106 << "\t-c - Curve" << std::endl
107 << "\t-V - Vertex" << std::endl
108 << "\t-m - Material set (block)" << std::endl
109 << "\t-d - Dirichlet set (nodeset)" << std::endl
110 << "\t-n - Neumann set (sideset)" << std::endl
111 << "\t-D - Parallel partitioning set (PARALLEL_PARTITION)" << std::endl
112 << "\tThe presence of one or more of the following flags limits " << std::endl
113 << "\tthe exported mesh to only elements of the corresponding " << std::endl
114 << "\tdimension. Vertices are always exported." << std::endl
115 << "\t-1 - Edges " << std::endl
116 << "\t-2 - Tri, Quad, Polygon " << std::endl
117 << "\t-3 - Tet, Hex, Prism, etc. " << std::endl;
118 }
119
120 static void print_help( const char* name )
121 {
122 std::cout << " This program can be used to convert between mesh file\n"
123 " formats, extract a subset of a mesh file to a separate\n"
124 " file, or both. The type of file to write is determined\n"
125 " from the file extension (e.g. \".vtk\") portion of the\n"
126 " output file name.\n"
127 " \n"
128 " While MOAB strives to export and import all data from\n"
129 " each supported file format, most file formats do\n"
130 " not support MOAB's entire data model. Thus MOAB cannot\n"
131 " guarantee lossless conversion for any file formats\n"
132 " other than the native HDF5 representation.\n"
133 "\n";
134
135 print_usage( name, std::cout );
136 exit( 0 );
137 }
138
139 static void usage_error( const char* name )
140 {
141 print_usage( name, std::cerr );
142 #ifdef MOAB_HAVE_MPI
143 MPI_Finalize();
144 #endif
145 exit( USAGE_ERROR );
146 }
147
148 static void list_formats( Interface* );
149 static bool parse_id_list( const char* string, std::set< int >& );
150 static void print_id_list( const char*, std::ostream& stream, const std::set< int >& list );
151 static void reset_times();
152 static void write_times( std::ostream& stream );
153 static void remove_entities_from_sets( Interface* gMB, Range& dead_entities, Range& empty_sets );
154 static void remove_from_vector( std::vector< EntityHandle >& vect, const Range& ents_to_remove );
155 static bool make_opts_string( std::vector< std::string > options, std::string& result );
156 static std::string percent_subst( const std::string& s, int val );
157
158 static int process_partition_file( Interface* gMB, std::string& metis_partition_file );
159
160 int main( int argc, char* argv[] )
161 {
162 int proc_id = 0;
163 #ifdef MOAB_HAVE_MPI
164 MPI_Init( &argc, &argv );
165 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
166 #endif
167
168 #ifdef MOAB_HAVE_TEMPESTREMAP
169 bool tempestin = false, tempestout = false;
170 #endif
171
172 Core core;
173 Interface* gMB = &core;
174 ErrorCode result;
175 Range range;
176
177 #if( defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_TEMPESTREMAP ) )
178 moab::ParallelComm* pcomm = new moab::ParallelComm( gMB, MPI_COMM_WORLD, 0 );
179 #endif
180
181 bool append_rank = false;
182 bool percent_rank_subst = false;
183 bool file_written = false;
184 int i, dim;
185 std::list< std::string >::iterator j;
186 bool dims[4] = { false, false, false, false };
187 const char* format = NULL; // output file format
188 std::list< std::string > in; // input file name list
189 std::string out; // output file name
190 bool verbose = false;
191 std::set< int > geom[4], mesh[4]; // user-specified IDs
192 std::vector< EntityHandle > set_list; // list of user-specified sets to write
193 std::vector< std::string > write_opts, read_opts;
194 std::string metis_partition_file;
195 #ifdef MOAB_HAVE_TEMPESTREMAP
196 std::string globalid_tag_name;
197 int spectral_order = 1;
198 bool unitscaling = false;
199 #endif
200
201 const char* const mesh_tag_names[] = { DIRICHLET_SET_TAG_NAME, NEUMANN_SET_TAG_NAME, MATERIAL_SET_TAG_NAME,
202 PARALLEL_PARTITION_TAG_NAME };
203 const char* const geom_names[] = { "VERTEX", "CURVE", "SURFACE", "VOLUME" };
204
205 // scan arguments
206 bool do_flag = true;
207 bool print_times = false;
208 bool generate[] = { false, false, false };
209 bool pval;
210 bool parallel = false, resolve_shared = false, exchange_ghosts = false;
211 for( i = 1; i < argc; i++ )
212 {
213 if( !argv[i][0] ) usage_error( argv[0] );
214
215 if( do_flag && argv[i][0] == '-' )
216 {
217 if( !argv[i][1] || ( argv[i][1] != 'M' && argv[i][2] ) ) usage_error( argv[0] );
218
219 switch( argv[i][1] )
220 {
221 // do flag arguments:
222 case '-':
223 do_flag = false;
224 break;
225 case 'g':
226 verbose = true;
227 break;
228 case 't':
229 print_times = true;
230 break;
231 case 'A':
232 break;
233 case 'h':
234 case 'H':
235 print_help( argv[0] );
236 break;
237 case 'l':
238 list_formats( gMB );
239 break;
240 #ifdef MOAB_HAVE_MPI
241 case 'P':
242 append_rank = true;
243 break;
244 case 'p':
245 percent_rank_subst = true;
246 break;
247 case 'M':
248 parallel = true;
249 if( argv[i][2] == '1' || argv[i][2] == '2' ) resolve_shared = true;
250 if( argv[i][2] == '2' ) exchange_ghosts = true;
251 break;
252 #endif
253 #ifdef MOAB_HAVE_TEMPESTREMAP
254 case 'B':
255 tempestin = true;
256 break;
257 case 'b':
258 tempestout = true;
259 break;
260 case 'S':
261 unitscaling = true;
262 break;
263 #endif
264 case '1':
265 case '2':
266 case '3':
267 dims[argv[i][1] - '0'] = true;
268 break;
269 // do options that require additional args:
270 default:
271 ++i;
272 if( i == argc || argv[i][0] == '-' )
273 {
274 std::cerr << "Expected argument following " << argv[i - 1] << std::endl;
275 usage_error( argv[0] );
276 }
277 if( argv[i - 1][1] == 'I' )
278 {
279 dim = atoi( argv[i] );
280 if( dim < 1 || dim > 2 )
281 {
282 std::cerr << "Invalid dimension value following -I" << std::endl;
283 usage_error( argv[0] );
284 }
285 generate[dim] = true;
286 continue;
287 }
288 pval = false;
289 switch( argv[i - 1][1] )
290 {
291 case 'a':
292 read_opts.push_back( std::string( "SAT_FILE=" ) + argv[i] );
293 pval = true;
294 break;
295 case 'f':
296 format = argv[i];
297 pval = true;
298 break;
299 case 'o':
300 write_opts.push_back( argv[i] );
301 pval = true;
302 break;
303 case 'O':
304 read_opts.push_back( argv[i] );
305 pval = true;
306 break;
307 #ifdef MOAB_HAVE_TEMPESTREMAP
308 case 'i':
309 globalid_tag_name = std::string( argv[i] );
310 pval = true;
311 break;
312 case 'r':
313 spectral_order = atoi( argv[i] );
314 pval = true;
315 break;
316 #endif
317 case 'v':
318 pval = parse_id_list( argv[i], geom[3] );
319 break;
320 case 's':
321 pval = parse_id_list( argv[i], geom[2] );
322 break;
323 case 'c':
324 pval = parse_id_list( argv[i], geom[1] );
325 break;
326 case 'V':
327 pval = parse_id_list( argv[i], geom[0] );
328 break;
329 case 'D':
330 pval = parse_id_list( argv[i], mesh[3] );
331 break;
332 case 'm':
333 pval = parse_id_list( argv[i], mesh[2] );
334 break;
335 case 'n':
336 pval = parse_id_list( argv[i], mesh[1] );
337 break;
338 case 'd':
339 pval = parse_id_list( argv[i], mesh[0] );
340 break;
341 case 'z':
342 metis_partition_file = argv[i];
343 pval = true;
344 break;
345 default:
346 std::cerr << "Invalid option: " << argv[i] << std::endl;
347 }
348
349 if( !pval )
350 {
351 std::cerr << "Invalid flag or flag value: " << argv[i - 1] << " " << argv[i] << std::endl;
352 usage_error( argv[0] );
353 }
354 }
355 }
356 // do file names
357 else
358 {
359 in.push_back( argv[i] );
360 }
361 }
362 if( in.size() < 2 )
363 {
364 std::cerr << "No output file name specified." << std::endl;
365 usage_error( argv[0] );
366 }
367 // output file name is the last one specified
368 out = in.back();
369 in.pop_back();
370
371 if( append_rank )
372 {
373 std::ostringstream mod;
374 mod << out << "." << proc_id;
375 out = mod.str();
376 }
377
378 if( percent_rank_subst )
379 {
380 for( j = in.begin(); j != in.end(); ++j )
381 *j = percent_subst( *j, proc_id );
382 out = percent_subst( out, proc_id );
383 }
384
385 // construct options string from individual options
386 std::string read_options, write_options;
387 if( parallel )
388 {
389 read_opts.push_back( "PARALLEL=READ_PART" );
390 read_opts.push_back( "PARTITION=PARALLEL_PARTITION" );
391 if( !append_rank && !percent_rank_subst ) write_opts.push_back( "PARALLEL=WRITE_PART" );
392 }
393 if( resolve_shared ) read_opts.push_back( "PARALLEL_RESOLVE_SHARED_ENTS" );
394 if( exchange_ghosts ) read_opts.push_back( "PARALLEL_GHOSTS=3.0.1" );
395
396 if( !make_opts_string( read_opts, read_options ) || !make_opts_string( write_opts, write_options ) )
397 {
398 #ifdef MOAB_HAVE_MPI
399 MPI_Finalize();
400 #endif
401 return USAGE_ERROR;
402 }
403
404 if( !metis_partition_file.empty() )
405 {
406 if( ( in.size() != 1 ) || ( proc_id != 0 ) )
407 {
408 std::cerr << " mpas partition allows only one input file, in serial conversion\n";
409 #ifdef MOAB_HAVE_MPI
410 MPI_Finalize();
411 #endif
412 return USAGE_ERROR;
413 }
414 }
415
416 Tag id_tag = gMB->globalId_tag();
417
418 // Read the input file.
419 #ifdef MOAB_HAVE_TEMPESTREMAP
420 if( tempestin && in.size() > 1 )
421 {
422 std::cerr << " we can read only one tempest files at a time\n";
423 #ifdef MOAB_HAVE_MPI
424 MPI_Finalize();
425 #endif
426 return USAGE_ERROR;
427 }
428 TempestRemapper* remapper = NULL;
429 if( tempestin or tempestout )
430 {
431 #ifdef MOAB_HAVE_MPI
432 remapper = new moab::TempestRemapper( gMB, pcomm, offlineGenerator );
433 #else
434 remapper = new moab::TempestRemapper( gMB, offlineGenerator );
435 #endif
436 }
437
438 bool use_overlap_context = false;
439 Tag srcParentTag, tgtParentTag;
440
441 #endif
442 for( j = in.begin(); j != in.end(); ++j )
443 {
444 std::string inFileName = *j;
445
446 reset_times();
447
448 #ifdef MOAB_HAVE_TEMPESTREMAP
449 if( remapper )
450 {
451 remapper->meshValidate = false;
452 // remapper->constructEdgeMap = true;
453 remapper->initialize();
454
455 if( tempestin )
456 {
457 moab::EntityHandle& srcmesh = remapper->GetMeshSet( moab::Remapper::SourceMesh );
458
459 // convert
460 result = remapper->LoadMesh( moab::Remapper::SourceMesh, inFileName, moab::TempestRemapper::DEFAULT );MB_CHK_ERR( result );
461
462 Mesh* tempestMesh = remapper->GetMesh( moab::Remapper::SourceMesh );
463 tempestMesh->RemoveZeroEdges();
464 tempestMesh->RemoveCoincidentNodes();
465
466 // Load the meshes and validate
467 result = remapper->ConvertTempestMesh( moab::Remapper::SourceMesh );
468
469 if( unitscaling )
470 {
471 result = moab::IntxUtils::ScaleToRadius( gMB, srcmesh, 1.0 );MB_CHK_ERR( result );
472 }
473
474 // Check if we are converting a RLL grid
475 NcFile ncInput( inFileName.c_str(), NcFile::ReadOnly );
476
477 NcError error_temp( NcError::silent_nonfatal );
478 // get the attribute
479 NcAtt* attRectilinear = ncInput.get_att( "rectilinear" );
480 NcVar* varGridDims = ncInput.get_var( "grid_dims" );
481
482 // If rectilinear attribute present, mark it
483 std::vector< int > vecDimSizes( 3, 0 );
484 Tag rectilinearTag;
485 // Tag data contains: guessed mesh type, mesh size1, mesh size 2
486 // Example: CS(0)/ICO(1)/ICOD(2), num_elements, num_nodes
487 // : RLL(3), num_lat, num_lon
488 result = gMB->tag_get_handle( "ClimateMetadata", 3, MB_TYPE_INTEGER, rectilinearTag,
489 MB_TAG_SPARSE | MB_TAG_CREAT, vecDimSizes.data() );MB_CHK_SET_ERR( result, "can't create rectilinear sizes tag" );
490
491 if( attRectilinear != nullptr )
492 {
493 // Obtain rectilinear attributes (dimension sizes)
494 NcAtt* attRectilinearDim0Size = ncInput.get_att( "rectilinear_dim0_size" );
495 NcAtt* attRectilinearDim1Size = ncInput.get_att( "rectilinear_dim1_size" );
496
497 if( attRectilinearDim0Size == nullptr )
498 {
499 _EXCEPTIONT( "Missing attribute \"rectilinear_dim0_size\"" );
500 }
501 if( attRectilinearDim1Size == nullptr )
502 {
503 _EXCEPTIONT( "Missing attribute \"rectilinear_dim1_size\"" );
504 }
505
506 int nDim0Size = attRectilinearDim0Size->as_int( 0 );
507 int nDim1Size = attRectilinearDim1Size->as_int( 0 );
508
509 // Obtain rectilinear attributes (dimension names)
510 NcAtt* attRectilinearDim0Name = ncInput.get_att( "rectilinear_dim0_name" );
511 NcAtt* attRectilinearDim1Name = ncInput.get_att( "rectilinear_dim1_name" );
512
513 if( attRectilinearDim0Name == nullptr )
514 {
515 _EXCEPTIONT( "Missing attribute \"rectilinear_dim0_name\"" );
516 }
517 if( attRectilinearDim1Name == nullptr )
518 {
519 _EXCEPTIONT( "Missing attribute \"rectilinear_dim1_name\"" );
520 }
521
522 std::string strDim0Name = attRectilinearDim0Name->as_string( 0 );
523 std::string strDim1Name = attRectilinearDim1Name->as_string( 0 );
524
525 std::map< std::string, int > vecDimNameSizes;
526 // Push rectilinear attributes into array
527 vecDimNameSizes[strDim0Name] = nDim0Size;
528 vecDimNameSizes[strDim1Name] = nDim1Size;
529 vecDimSizes[0] = static_cast< int >( moab::TempestRemapper::RLL );
530 vecDimSizes[1] = vecDimNameSizes["lat"];
531 vecDimSizes[2] = vecDimNameSizes["lon"];
532 }
533 else if( varGridDims != nullptr )
534 {
535 // Obtain rectilinear attributes (dimension sizes)
536 NcDim* dimGridRank = varGridDims->get_dim( 0 );
537 if( dimGridRank == NULL )
538 {
539 _EXCEPTIONT( "Variable \"grid_dims\" has no dimensions" );
540 }
541
542 int gridrank = dimGridRank->size();
543 // now get the rank dimensions
544 int gridsizes[2];
545 varGridDims->get( &( gridsizes[0] ), dimGridRank->size() );
546
547 const Range& elems = remapper->GetMeshEntities( moab::Remapper::SourceMesh );
548 bool isQuads = elems.all_of_type( moab::MBQUAD );
549 bool isTris = elems.all_of_type( moab::MBTRI );
550
551 if( gridrank == 1 )
552 {
553 vecDimSizes[0] = ( isQuads ? static_cast< int >( moab::TempestRemapper::CS )
554 : ( isTris ? static_cast< int >( moab::TempestRemapper::ICO )
555 : static_cast< int >( moab::TempestRemapper::ICOD ) ) );
556 vecDimSizes[1] = elems.size();
557 vecDimSizes[2] = remapper->GetMeshVertices( moab::Remapper::SourceMesh ).size();
558 }
559 else
560 {
561 vecDimSizes[0] = static_cast< int >( moab::TempestRemapper::RLL );
562 vecDimSizes[1] = gridsizes[0];
563 vecDimSizes[2] = gridsizes[1];
564 }
565 }
566 else
567 {
568 const Range& elems = remapper->GetMeshEntities( moab::Remapper::SourceMesh );
569 bool isQuads = elems.all_of_type( moab::MBQUAD );
570 bool isTris = elems.all_of_type( moab::MBTRI );
571 // vecDimSizes[0] = static_cast< int >( remapper->GetMeshType( moab::Remapper::SourceMesh );
572 vecDimSizes[0] = ( isQuads ? static_cast< int >( moab::TempestRemapper::CS )
573 : ( isTris ? static_cast< int >( moab::TempestRemapper::ICO )
574 : static_cast< int >( moab::TempestRemapper::ICOD ) ) );
575 vecDimSizes[1] = elems.size();
576 vecDimSizes[2] = remapper->GetMeshVertices( moab::Remapper::SourceMesh ).size();
577 }
578 moab::EntityHandle mSet = 0;
579 // mSet = remapper->GetMeshSet( moab::Remapper::SourceMesh );
580 result = gMB->tag_set_data( rectilinearTag, &mSet, 1, vecDimSizes.data() );MB_CHK_ERR( result );
581
582 switch( vecDimSizes[0] )
583 {
584 case 0:
585 printf( "Cubed-Sphere mesh: %d (elements), %d (vertices)\n", vecDimSizes[1], vecDimSizes[2] );
586 break;
587 case 1:
588 printf( "Rectilinear RLL mesh: (lon) %d X (lat) %d\n", vecDimSizes[2], vecDimSizes[1] );
589 break;
590 case 2:
591 printf( "Icosahedral (triangular) mesh: %d (elements), %d (vertices)\n", vecDimSizes[1],
592 vecDimSizes[2] );
593 break;
594 case 3:
595 default:
596 printf( "Polygonal mesh: %d (elements), %d (vertices)\n", vecDimSizes[1], vecDimSizes[2] );
597 break;
598 }
599
600 ncInput.close();
601
602 const size_t nOverlapFaces = tempestMesh->faces.size();
603 if( tempestMesh->vecSourceFaceIx.size() == nOverlapFaces &&
604 tempestMesh->vecSourceFaceIx.size() == nOverlapFaces )
605 {
606 int defaultInt = -1;
607 use_overlap_context = true;
608 // Check if our MOAB mesh has RED and BLUE tags; this would indicate we are
609 // converting an overlap grid
610 result = gMB->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag,
611 MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( result, "can't create target parent tag" );
612
613 result = gMB->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag,
614 MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( result, "can't create source parent tag" );
615
616 const Range& faces = remapper->GetMeshEntities( moab::Remapper::SourceMesh );
617
618 std::vector< int > gids( faces.size() ), srcpar( faces.size() ), tgtpar( faces.size() );
619 result = gMB->tag_get_data( id_tag, faces, &gids[0] );MB_CHK_ERR( result );
620
621 for( unsigned ii = 0; ii < faces.size(); ++ii )
622 {
623 srcpar[ii] = tempestMesh->vecSourceFaceIx[gids[ii] - 1];
624 tgtpar[ii] = tempestMesh->vecTargetFaceIx[gids[ii] - 1];
625 }
626
627 result = gMB->tag_set_data( srcParentTag, faces, &srcpar[0] );MB_CHK_ERR( result );
628 result = gMB->tag_set_data( tgtParentTag, faces, &tgtpar[0] );MB_CHK_ERR( result );
629
630 srcpar.clear();
631 tgtpar.clear();
632 gids.clear();
633 }
634 }
635 else if( tempestout )
636 {
637 moab::EntityHandle& srcmesh = remapper->GetMeshSet( moab::Remapper::SourceMesh );
638 moab::EntityHandle& ovmesh = remapper->GetMeshSet( moab::Remapper::OverlapMesh );
639
640 // load the mesh in MOAB format
641 std::vector< int > metadata( 2 );
642 result = remapper->LoadNativeMesh( *j, srcmesh, metadata );MB_CHK_ERR( result );
643
644 if( unitscaling )
645 {
646 result = moab::IntxUtils::ScaleToRadius( gMB, srcmesh, 1.0 );MB_CHK_ERR( result );
647 }
648
649 // Check if our MOAB mesh has RED and BLUE tags; this would indicate we are converting
650 // an overlap grid
651 ErrorCode rval1 = gMB->tag_get_handle( "SourceParent", srcParentTag );
652 ErrorCode rval2 = gMB->tag_get_handle( "TargetParent", tgtParentTag );
653 if( rval1 == MB_SUCCESS && rval2 == MB_SUCCESS )
654 {
655 use_overlap_context = true;
656 ovmesh = srcmesh;
657
658 Tag countTag;
659 result = gMB->tag_get_handle( "Counting", countTag );MB_CHK_ERR( result );
660
661 // Load the meshes and validate
662 Tag order;
663 ReorderTool reorder_tool( &core );
664 result = reorder_tool.handle_order_from_int_tag( srcParentTag, -1, order );MB_CHK_ERR( result );
665 result = reorder_tool.reorder_entities( order );MB_CHK_ERR( result );
666 result = gMB->tag_delete( order );MB_CHK_ERR( result );
667 result = remapper->ConvertMeshToTempest( moab::Remapper::OverlapMesh );MB_CHK_ERR( result );
668 }
669 else
670 {
671 if( metadata[0] == static_cast< int >( moab::TempestRemapper::RLL ) )
672 {
673 assert( metadata.size() );
674 std::cout << "Converting a RLL mesh with rectilinear dimension: " << metadata[0] << " X "
675 << metadata[1] << std::endl;
676 }
677
678 // Convert the mesh and validate
679 result = remapper->ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( result );
680 }
681 }
682 }
683 else
684 result = gMB->load_file( j->c_str(), 0, read_options.c_str() );
685 #else
686 result = gMB->load_file( j->c_str(), 0, read_options.c_str() );
687 #endif
688 if( MB_SUCCESS != result )
689 {
690 std::cerr << "Failed to load \"" << *j << "\"." << std::endl;
691 std::cerr << "Error code: " << gMB->get_error_string( result ) << " (" << result << ")" << std::endl;
692 std::string message;
693 if( MB_SUCCESS == gMB->get_last_error( message ) && !message.empty() )
694 std::cerr << "Error message: " << message << std::endl;
695 #ifdef MOAB_HAVE_MPI
696 MPI_Finalize();
697 #endif
698 return READ_ERROR;
699 }
700 if( !proc_id ) std::cerr << "Read \"" << *j << "\"" << std::endl;
701 if( print_times && !proc_id ) write_times( std::cout );
702 }
703
704 // Determine if the user has specified any geometry sets to write
705 bool have_geom = false;
706 for( dim = 0; dim <= 3; ++dim )
707 {
708 if( !geom[dim].empty() ) have_geom = true;
709 if( verbose ) print_id_list( geom_names[dim], std::cout, geom[dim] );
710 }
711
712 // True if the user has specified any sets to write
713 bool have_sets = have_geom;
714
715 // Get geometry tags
716 Tag dim_tag;
717 if( have_geom )
718 {
719 if( id_tag == 0 )
720 {
721 std::cerr << "No ID tag defined." << std::endl;
722 have_geom = false;
723 }
724 result = gMB->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dim_tag );
725 if( MB_SUCCESS != result )
726 {
727 std::cerr << "No geometry tag defined." << std::endl;
728 have_geom = false;
729 }
730 }
731
732 // Get geometry sets
733 if( have_geom )
734 {
735 int id_val;
736 Tag tags[] = { id_tag, dim_tag };
737 const void* vals[] = { &id_val, &dim };
738 for( dim = 0; dim <= 3; ++dim )
739 {
740 int init_count = set_list.size();
741 for( std::set< int >::iterator iter = geom[dim].begin(); iter != geom[dim].end(); ++iter )
742 {
743 id_val = *iter;
744 range.clear();
745 result = gMB->get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, range );
746 if( MB_SUCCESS != result || range.empty() )
747 {
748 range.clear();
749 std::cerr << geom_names[dim] << " " << id_val << " not found.\n";
750 }
751 std::copy( range.begin(), range.end(), std::back_inserter( set_list ) );
752 }
753
754 if( verbose )
755 std::cout << "Found " << ( set_list.size() - init_count ) << ' ' << geom_names[dim] << " sets"
756 << std::endl;
757 }
758 }
759
760 // Get mesh groupings
761 for( i = 0; i < 4; ++i )
762 {
763 if( verbose ) print_id_list( mesh_tag_names[i], std::cout, mesh[i] );
764
765 if( mesh[i].empty() ) continue;
766 have_sets = true;
767
768 // Get tag
769 Tag tag;
770 result = gMB->tag_get_handle( mesh_tag_names[i], 1, MB_TYPE_INTEGER, tag );
771 if( MB_SUCCESS != result )
772 {
773 std::cerr << "Tag not found: " << mesh_tag_names[i] << std::endl;
774 continue;
775 }
776
777 // get entity sets
778 int init_count = set_list.size();
779 for( std::set< int >::iterator iter = mesh[i].begin(); iter != mesh[i].end(); ++iter )
780 {
781 range.clear();
782 const void* vals[] = { &*iter };
783 result = gMB->get_entities_by_type_and_tag( 0, MBENTITYSET, &tag, vals, 1, range );
784 if( MB_SUCCESS != result || range.empty() )
785 {
786 range.clear();
787 std::cerr << mesh_tag_names[i] << " " << *iter << " not found.\n";
788 }
789 std::copy( range.begin(), range.end(), std::back_inserter( set_list ) );
790 }
791
792 if( verbose )
793 std::cout << "Found " << ( set_list.size() - init_count ) << ' ' << mesh_tag_names[i] << " sets"
794 << std::endl;
795 }
796
797 // Check if output is limited to certain dimensions of elements
798 bool bydim = false;
799 for( dim = 1; dim < 4; ++dim )
800 if( dims[dim] ) bydim = true;
801
802 // Check conflicting input
803 if( bydim )
804 {
805 if( generate[1] && !dims[1] )
806 {
807 std::cerr << "Warning: Request to generate 1D internal entities but not export them." << std::endl;
808 generate[1] = false;
809 }
810 if( generate[2] && !dims[2] )
811 {
812 std::cerr << "Warning: Request to generate 2D internal entities but not export them." << std::endl;
813 generate[2] = false;
814 }
815 }
816
817 // Generate any internal entities
818 if( generate[1] || generate[2] )
819 {
820 EntityHandle all_mesh = 0;
821 const EntityHandle* sets = &all_mesh;
822 int num_sets = 1;
823 if( have_sets )
824 {
825 num_sets = set_list.size();
826 sets = &set_list[0];
827 }
828 for( i = 0; i < num_sets; ++i )
829 {
830 Range dim3, dim2, adj;
831 gMB->get_entities_by_dimension( sets[i], 3, dim3, true );
832 if( generate[1] )
833 {
834 gMB->get_entities_by_dimension( sets[i], 2, dim2, true );
835 gMB->get_adjacencies( dim3, 1, true, adj, Interface::UNION );
836 gMB->get_adjacencies( dim2, 1, true, adj, Interface::UNION );
837 }
838 if( generate[2] )
839 {
840 gMB->get_adjacencies( dim3, 2, true, adj, Interface::UNION );
841 }
842 if( sets[i] ) gMB->add_entities( sets[i], adj );
843 }
844 }
845
846 // Delete any entities not of the dimensions to be exported
847 if( bydim )
848 {
849 // Get list of dead elements
850 Range dead_entities, tmp_range;
851 for( dim = 1; dim <= 3; ++dim )
852 {
853 if( dims[dim] ) continue;
854 gMB->get_entities_by_dimension( 0, dim, tmp_range );
855 dead_entities.merge( tmp_range );
856 }
857 // Remove dead entities from all sets, and add all
858 // empty sets to list of dead entities.
859 Range empty_sets;
860 remove_entities_from_sets( gMB, dead_entities, empty_sets );
861 while( !empty_sets.empty() )
862 {
863 if( !set_list.empty() ) remove_from_vector( set_list, empty_sets );
864 dead_entities.merge( empty_sets );
865 tmp_range.clear();
866 remove_entities_from_sets( gMB, empty_sets, tmp_range );
867 empty_sets = subtract( tmp_range, dead_entities );
868 }
869 // Destroy dead entities
870 gMB->delete_entities( dead_entities );
871 }
872
873 // If user specified sets to write, but none were found, exit.
874 if( have_sets && set_list.empty() )
875 {
876 std::cerr << "Nothing to write." << std::endl;
877 #ifdef MOAB_HAVE_MPI
878 MPI_Finalize();
879 #endif
880 return ENT_NOT_FOUND;
881 }
882
883 // interpret the mpas partition file created by gpmetis
884 if( !metis_partition_file.empty() )
885 {
886 int err = process_partition_file( gMB, metis_partition_file );
887 if( err )
888 {
889 std::cerr << "Failed to process partition file \"" << metis_partition_file << "\"." << std::endl;
890 #ifdef MOAB_HAVE_MPI
891 MPI_Finalize();
892 #endif
893 return WRITE_ERROR;
894 }
895 }
896 if( verbose )
897 {
898 if( have_sets )
899 std::cout << "Found " << set_list.size() << " specified sets to write (total)." << std::endl;
900 else
901 std::cout << "No sets specifed. Writing entire mesh." << std::endl;
902 }
903
904 // Write the output file
905 reset_times();
906 #ifdef MOAB_HAVE_TEMPESTREMAP
907 if( remapper )
908 {
909 Range faces;
910 Mesh* tempestMesh =
911 remapper->GetMesh( ( use_overlap_context ? moab::Remapper::OverlapMesh : moab::Remapper::SourceMesh ) );
912 moab::EntityHandle& srcmesh =
913 remapper->GetMeshSet( ( use_overlap_context ? moab::Remapper::OverlapMesh : moab::Remapper::SourceMesh ) );
914 result = gMB->get_entities_by_dimension( srcmesh, 2, faces );MB_CHK_ERR( result );
915 int ntot_elements = 0, nelements = faces.size();
916 #ifdef MOAB_HAVE_MPI
917 int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, pcomm->comm() );
918 if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
919 #else
920 ntot_elements = nelements;
921 #endif
922
923 Tag gidTag = gMB->globalId_tag();
924 std::vector< int > gids( faces.size() );
925 result = gMB->tag_get_data( gidTag, faces, &gids[0] );MB_CHK_ERR( result );
926
927 if( faces.size() > 1 && gids[0] == gids[1] && !use_overlap_context )
928 {
929 #ifdef MOAB_HAVE_MPI
930 result = pcomm->assign_global_ids( srcmesh, 2, 1, false );MB_CHK_ERR( result );
931 #else
932 result = remapper->assign_vertex_element_IDs( gidTag, srcmesh, 2, 1 );MB_CHK_ERR( result );
933 result = remapper->assign_vertex_element_IDs( gidTag, srcmesh, 0, 1 );MB_CHK_ERR( result );
934 #endif
935 }
936
937 // VSM: If user requested explicitly for some metadata, we need to generate the DoF ID tag
938 // and set the appropriate numbering based on specified discretization order
939 // Useful only for SE meshes with GLL DoFs
940 if( spectral_order > 1 && globalid_tag_name.size() > 1 )
941 {
942 result = remapper->GenerateMeshMetadata( *tempestMesh, ntot_elements, faces, NULL, globalid_tag_name,
943 spectral_order );MB_CHK_ERR( result );
944 }
945
946 if( tempestout )
947 {
948 // Check if our MOAB mesh has RED and BLUE tags; this would indicate we are converting an
949 // overlap grid
950 if( use_overlap_context && false )
951 {
952 const int nOverlapFaces = faces.size();
953 // Overlap mesh: resize the source and target connection arrays
954 tempestMesh->vecSourceFaceIx.resize( nOverlapFaces ); // 0-based indices corresponding to source mesh
955 tempestMesh->vecTargetFaceIx.resize( nOverlapFaces ); // 0-based indices corresponding to target mesh
956 result = gMB->tag_get_data( srcParentTag, faces, &tempestMesh->vecSourceFaceIx[0] );MB_CHK_ERR( result );
957 result = gMB->tag_get_data( tgtParentTag, faces, &tempestMesh->vecTargetFaceIx[0] );MB_CHK_ERR( result );
958 }
959 // Write out the mesh using TempestRemap
960 tempestMesh->Write( out, NcFile::Netcdf4 );
961 file_written = true;
962 }
963 delete remapper; // cleanup
964 }
965 #endif
966
967 if( !file_written )
968 {
969 if( have_sets )
970 result = gMB->write_file( out.c_str(), format, write_options.c_str(), &set_list[0], set_list.size() );
971 else
972 result = gMB->write_file( out.c_str(), format, write_options.c_str() );
973 if( MB_SUCCESS != result )
974 {
975 std::cerr << "Failed to write \"" << out << "\"." << std::endl;
976 std::cerr << "Error code: " << gMB->get_error_string( result ) << " (" << result << ")" << std::endl;
977 std::string message;
978 if( MB_SUCCESS == gMB->get_last_error( message ) && !message.empty() )
979 std::cerr << "Error message: " << message << std::endl;
980 #ifdef MOAB_HAVE_MPI
981 MPI_Finalize();
982 #endif
983 return WRITE_ERROR;
984 }
985 }
986
987 if( !proc_id ) std::cerr << "Wrote \"" << out << "\"" << std::endl;
988 if( print_times && !proc_id ) write_times( std::cout );
989
990 #ifdef MOAB_HAVE_MPI
991 MPI_Finalize();
992 #endif
993 return 0;
994 }
995
996 bool parse_id_list( const char* string, std::set< int >& results )
997 {
998 bool okay = true;
999 char* mystr = strdup( string );
1000 for( const char* ptr = strtok( mystr, "," ); ptr; ptr = strtok( 0, "," ) )
1001 {
1002 char* endptr;
1003 long val = strtol( ptr, &endptr, 0 );
1004 if( endptr == ptr || val <= 0 )
1005 {
1006 std::cerr << "Not a valid id: " << ptr << std::endl;
1007 okay = false;
1008 break;
1009 }
1010
1011 long val2 = val;
1012 if( *endptr == '-' )
1013 {
1014 const char* sptr = endptr + 1;
1015 val2 = strtol( sptr, &endptr, 0 );
1016 if( endptr == sptr || val2 <= 0 )
1017 {
1018 std::cerr << "Not a valid id: " << sptr << std::endl;
1019 okay = false;
1020 break;
1021 }
1022 if( val2 < val )
1023 {
1024 std::cerr << "Invalid id range: " << ptr << std::endl;
1025 okay = false;
1026 break;
1027 }
1028 }
1029
1030 if( *endptr )
1031 {
1032 std::cerr << "Unexpected character: " << *endptr << std::endl;
1033 okay = false;
1034 break;
1035 }
1036
1037 for( ; val <= val2; ++val )
1038 if( !results.insert( (int)val ).second ) std::cerr << "Warning: duplicate Id: " << val << std::endl;
1039 }
1040
1041 free( mystr );
1042 return okay;
1043 }
1044
1045 void print_id_list( const char* head, std::ostream& stream, const std::set< int >& list )
1046 {
1047 stream << head << ": ";
1048
1049 if( list.empty() )
1050 {
1051 stream << "(none)" << std::endl;
1052 return;
1053 }
1054
1055 int start, prev;
1056 std::set< int >::const_iterator iter = list.begin();
1057 start = prev = *( iter++ );
1058 for( ;; )
1059 {
1060 if( iter == list.end() || *iter != 1 + prev )
1061 {
1062 stream << start;
1063 if( prev != start ) stream << '-' << prev;
1064 if( iter == list.end() ) break;
1065 stream << ", ";
1066 start = *iter;
1067 }
1068 prev = *( iter++ );
1069 }
1070
1071 stream << std::endl;
1072 }
1073
1074 static void print_time( int clk_per_sec, const char* prefix, clock_t ticks, std::ostream& stream )
1075 {
1076 ticks *= clk_per_sec / 100;
1077 clock_t centi = ticks % 100;
1078 clock_t seconds = ticks / 100;
1079 stream << prefix;
1080 if( seconds < 120 )
1081 {
1082 stream << ( ticks / 100 ) << "." << centi << "s" << std::endl;
1083 }
1084 else
1085 {
1086 clock_t minutes = ( seconds / 60 ) % 60;
1087 clock_t hours = ( seconds / 3600 );
1088 seconds %= 60;
1089 if( hours ) stream << hours << "h";
1090 if( minutes ) stream << minutes << "m";
1091 if( seconds || centi ) stream << seconds << "." << centi << "s";
1092 stream << " (" << ( ticks / 100 ) << "." << centi << "s)" << std::endl;
1093 }
1094 }
1095
1096 clock_t usr_time, sys_time, abs_time;
1097
1098 #ifdef WIN32
1099
1100 void reset_times()
1101 {
1102 abs_time = clock();
1103 }
1104
1105 void write_times( std::ostream& stream )
1106 {
1107 clock_t abs_tm = clock();
1108 print_time( CLOCKS_PER_SEC, " ", abs_tm - abs_time, stream );
1109 abs_time = abs_tm;
1110 }
1111
1112 #else
1113
1114 void reset_times()
1115 {
1116 tms timebuf;
1117 abs_time = times( &timebuf );
1118 usr_time = timebuf.tms_utime;
1119 sys_time = timebuf.tms_stime;
1120 }
1121
1122 void write_times( std::ostream& stream )
1123 {
1124 clock_t usr_tm, sys_tm, abs_tm;
1125 tms timebuf;
1126 abs_tm = times( &timebuf );
1127 usr_tm = timebuf.tms_utime;
1128 sys_tm = timebuf.tms_stime;
1129 print_time( sysconf( _SC_CLK_TCK ), " real: ", abs_tm - abs_time, stream );
1130 print_time( sysconf( _SC_CLK_TCK ), " user: ", usr_tm - usr_time, stream );
1131 print_time( sysconf( _SC_CLK_TCK ), " system: ", sys_tm - sys_time, stream );
1132 abs_time = abs_tm;
1133 usr_time = usr_tm;
1134 sys_time = sys_tm;
1135 }
1136
1137 #endif
1138
1139 bool make_opts_string( std::vector< std::string > options, std::string& opts )
1140 {
1141 opts.clear();
1142 if( options.empty() ) return true;
1143
1144 // choose a separator character
1145 std::vector< std::string >::const_iterator i;
1146 char separator = '\0';
1147 const char* alt_separators = ";+,:\t\n";
1148 for( const char* sep_ptr = alt_separators; *sep_ptr; ++sep_ptr )
1149 {
1150 bool seen = false;
1151 for( i = options.begin(); i != options.end(); ++i )
1152 if( i->find( *sep_ptr, 0 ) != std::string::npos )
1153 {
1154 seen = true;
1155 break;
1156 }
1157 if( !seen )
1158 {
1159 separator = *sep_ptr;
1160 break;
1161 }
1162 }
1163 if( !separator )
1164 {
1165 std::cerr << "Error: cannot find separator character for options string" << std::endl;
1166 return false;
1167 }
1168 if( separator != ';' )
1169 {
1170 opts = ";";
1171 opts += separator;
1172 }
1173
1174 // concatenate options
1175 i = options.begin();
1176 opts += *i;
1177 for( ++i; i != options.end(); ++i )
1178 {
1179 opts += separator;
1180 opts += *i;
1181 }
1182
1183 return true;
1184 }
1185
1186 void list_formats( Interface* gMB )
1187 {
1188 const char iface_name[] = "ReaderWriterSet";
1189 ErrorCode err;
1190 ReaderWriterSet* set = 0;
1191 ReaderWriterSet::iterator i;
1192 std::ostream& str = std::cout;
1193
1194 // get ReaderWriterSet
1195 err = gMB->query_interface( set );
1196 if( err != MB_SUCCESS || !set )
1197 {
1198 std::cerr << "Internal error: Interface \"" << iface_name << "\" not available.\n";
1199 exit( OTHER_ERROR );
1200 }
1201
1202 // get field width for format description
1203 size_t w = 0;
1204 for( i = set->begin(); i != set->end(); ++i )
1205 if( i->description().length() > w ) w = i->description().length();
1206
1207 // write table header
1208 str << "Format " << std::setw( w ) << std::left << "Description" << " Read Write File Name Suffixes\n"
1209 << "------ " << std::setw( w ) << std::setfill( '-' ) << "" << std::setfill( ' ' )
1210 << " ---- ----- ------------------\n";
1211
1212 // write table data
1213 for( i = set->begin(); i != set->end(); ++i )
1214 {
1215 std::vector< std::string > ext;
1216 i->get_extensions( ext );
1217 str << std::setw( 6 ) << i->name() << " " << std::setw( w ) << std::left << i->description() << " "
1218 << ( i->have_reader() ? " yes" : " no" ) << " " << ( i->have_writer() ? " yes" : " no" ) << " ";
1219 for( std::vector< std::string >::iterator j = ext.begin(); j != ext.end(); ++j )
1220 str << " " << *j;
1221 str << std::endl;
1222 }
1223 str << std::endl;
1224
1225 gMB->release_interface( set );
1226 exit( 0 );
1227 }
1228
1229 void remove_entities_from_sets( Interface* gMB, Range& dead_entities, Range& empty_sets )
1230 {
1231 empty_sets.clear();
1232 Range sets;
1233 gMB->get_entities_by_type( 0, MBENTITYSET, sets );
1234 for( Range::iterator i = sets.begin(); i != sets.end(); ++i )
1235 {
1236 Range set_contents;
1237 gMB->get_entities_by_handle( *i, set_contents, false );
1238 set_contents = intersect( set_contents, dead_entities );
1239 gMB->remove_entities( *i, set_contents );
1240 set_contents.clear();
1241 gMB->get_entities_by_handle( *i, set_contents, false );
1242 if( set_contents.empty() ) empty_sets.insert( *i );
1243 }
1244 }
1245
1246 void remove_from_vector( std::vector< EntityHandle >& vect, const Range& ents_to_remove )
1247 {
1248 Range::const_iterator i;
1249 std::vector< EntityHandle >::iterator j;
1250 for( i = ents_to_remove.begin(); i != ents_to_remove.end(); ++i )
1251 {
1252 j = std::find( vect.begin(), vect.end(), *i );
1253 if( j != vect.end() ) vect.erase( j );
1254 }
1255 }
1256
1257 std::string percent_subst( const std::string& s, int val )
1258 {
1259 if( s.empty() ) return s;
1260
1261 size_t j = s.find( '%' );
1262 if( j == std::string::npos ) return s;
1263
1264 std::ostringstream st;
1265 st << s.substr( 0, j );
1266 st << val;
1267
1268 size_t i;
1269 while( ( i = s.find( '%', j + 1 ) ) != std::string::npos )
1270 {
1271 st << s.substr( j, i - j );
1272 st << val;
1273 j = i;
1274 }
1275 st << s.substr( j + 1 );
1276 return st.str();
1277 }
1278
1279 int process_partition_file( Interface* mb, std::string& metis_partition_file )
1280 {
1281 // how many faces in the file ? how do we make sure it is an mpas file?
1282 // mpas atmosphere files can be downloaded from here
1283 // https://mpas-dev.github.io/atmosphere/atmosphere_meshes.html
1284 Range faces;
1285 ErrorCode rval = mb->get_entities_by_dimension( 0, 2, faces );MB_CHK_ERR( rval );
1286 std::cout << " MPAS model has " << faces.size() << " polygons\n";
1287
1288 // read the partition file
1289 std::ifstream partfile;
1290 partfile.open( metis_partition_file.c_str() );
1291 std::string line;
1292 std::vector< int > parts;
1293 parts.resize( faces.size(), -1 );
1294 int i = 0;
1295 if( partfile.is_open() )
1296 {
1297 while( getline( partfile, line ) )
1298 {
1299 // cout << line << '\n';
1300 parts[i++] = atoi( line.c_str() );
1301 if( i > (int)faces.size() )
1302 {
1303 std::cerr << " too many lines in partition file \n. bail out \n";
1304 return 1;
1305 }
1306 }
1307 partfile.close();
1308 }
1309 std::vector< int >::iterator pmax = max_element( parts.begin(), parts.end() );
1310 std::vector< int >::iterator pmin = min_element( parts.begin(), parts.end() );
1311 if( *pmin <= -1 )
1312 {
1313 std::cerr << " partition file is incomplete, *pmin is -1 !! \n";
1314 return 1;
1315 }
1316 std::cout << " partitions range: " << *pmin << " " << *pmax << "\n";
1317 Tag part_set_tag;
1318 int dum_id = -1;
1319 rval = mb->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_set_tag, MB_TAG_SPARSE | MB_TAG_CREAT,
1320 &dum_id );MB_CHK_ERR( rval );
1321
1322 // get any sets already with this tag, and clear them
1323 // remove the parallel partition sets if they exist
1324 Range tagged_sets;
1325 rval = mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &part_set_tag, NULL, 1, tagged_sets, Interface::UNION );MB_CHK_ERR( rval );
1326 if( !tagged_sets.empty() )
1327 {
1328 rval = mb->clear_meshset( tagged_sets );MB_CHK_ERR( rval );
1329 rval = mb->tag_delete_data( part_set_tag, tagged_sets );MB_CHK_ERR( rval );
1330 }
1331 Tag gid;
1332 rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_ERR( rval );
1333 int num_sets = *pmax + 1;
1334 if( *pmin != 0 )
1335 {
1336 std::cout << " problem reading parts; min is not 0 \n";
1337 return 1;
1338 }
1339 for( i = 0; i < num_sets; i++ )
1340 {
1341 EntityHandle new_set;
1342 rval = mb->create_meshset( MESHSET_SET, new_set );MB_CHK_ERR( rval );
1343 tagged_sets.insert( new_set );
1344 }
1345 int* dum_ids = new int[num_sets];
1346 for( i = 0; i < num_sets; i++ )
1347 dum_ids[i] = i;
1348
1349 rval = mb->tag_set_data( part_set_tag, tagged_sets, dum_ids );MB_CHK_ERR( rval );
1350 delete[] dum_ids;
1351
1352 std::vector< int > gids;
1353 int num_faces = (int)faces.size();
1354 gids.resize( num_faces );
1355 rval = mb->tag_get_data( gid, faces, &gids[0] );MB_CHK_ERR( rval );
1356
1357 for( int j = 0; j < num_faces; j++ )
1358 {
1359 int eid = gids[j];
1360 EntityHandle eh = faces[j];
1361 int partition = parts[eid - 1];
1362 if( partition < 0 || partition >= num_sets )
1363 {
1364 std::cout << " wrong partition number \n";
1365 return 1;
1366 }
1367 rval = mb->add_entities( tagged_sets[partition], &eh, 1 );MB_CHK_ERR( rval );
1368 }
1369 return 0;
1370 }
1.9.1.