Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
convert.cpp
Go to the documentation of this file.
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 );
461  MB_CHK_ERR( result );
462 
463  Mesh* tempestMesh = remapper->GetMesh( moab::Remapper::SourceMesh );
464  tempestMesh->RemoveZeroEdges();
465  tempestMesh->RemoveCoincidentNodes();
466 
467  // Load the meshes and validate
468  result = remapper->ConvertTempestMesh( moab::Remapper::SourceMesh );
469 
470  if( unitscaling )
471  {
472  result = moab::IntxUtils::ScaleToRadius( gMB, srcmesh, 1.0 );
473  MB_CHK_ERR( result );
474  }
475 
476  // Check if we are converting a RLL grid
477  NcFile ncInput( inFileName.c_str(), NcFile::ReadOnly );
478 
479  NcError error_temp( NcError::silent_nonfatal );
480  // get the attribute
481  NcAtt* attRectilinear = ncInput.get_att( "rectilinear" );
482  NcVar* varGridDims = ncInput.get_var( "grid_dims" );
483 
484  // If rectilinear attribute present, mark it
485  std::vector< int > vecDimSizes( 3, 0 );
486  Tag rectilinearTag;
487  // Tag data contains: guessed mesh type, mesh size1, mesh size 2
488  // Example: CS(0)/ICO(1)/ICOD(2), num_elements, num_nodes
489  // : RLL(3), num_lat, num_lon
490  result = gMB->tag_get_handle( "ClimateMetadata", 3, MB_TYPE_INTEGER, rectilinearTag,
491  MB_TAG_SPARSE | MB_TAG_CREAT, vecDimSizes.data() );
492  MB_CHK_SET_ERR( result, "can't create rectilinear sizes tag" );
493 
494  if( attRectilinear != nullptr )
495  {
496  // Obtain rectilinear attributes (dimension sizes)
497  NcAtt* attRectilinearDim0Size = ncInput.get_att( "rectilinear_dim0_size" );
498  NcAtt* attRectilinearDim1Size = ncInput.get_att( "rectilinear_dim1_size" );
499 
500  if( attRectilinearDim0Size == nullptr )
501  {
502  _EXCEPTIONT( "Missing attribute \"rectilinear_dim0_size\"" );
503  }
504  if( attRectilinearDim1Size == nullptr )
505  {
506  _EXCEPTIONT( "Missing attribute \"rectilinear_dim1_size\"" );
507  }
508 
509  int nDim0Size = attRectilinearDim0Size->as_int( 0 );
510  int nDim1Size = attRectilinearDim1Size->as_int( 0 );
511 
512  // Obtain rectilinear attributes (dimension names)
513  NcAtt* attRectilinearDim0Name = ncInput.get_att( "rectilinear_dim0_name" );
514  NcAtt* attRectilinearDim1Name = ncInput.get_att( "rectilinear_dim1_name" );
515 
516  if( attRectilinearDim0Name == nullptr )
517  {
518  _EXCEPTIONT( "Missing attribute \"rectilinear_dim0_name\"" );
519  }
520  if( attRectilinearDim1Name == nullptr )
521  {
522  _EXCEPTIONT( "Missing attribute \"rectilinear_dim1_name\"" );
523  }
524 
525  std::string strDim0Name = attRectilinearDim0Name->as_string( 0 );
526  std::string strDim1Name = attRectilinearDim1Name->as_string( 0 );
527 
528  std::map< std::string, int > vecDimNameSizes;
529  // Push rectilinear attributes into array
530  vecDimNameSizes[strDim0Name] = nDim0Size;
531  vecDimNameSizes[strDim1Name] = nDim1Size;
532  vecDimSizes[0] = static_cast< int >( moab::TempestRemapper::RLL );
533  vecDimSizes[1] = vecDimNameSizes["lat"];
534  vecDimSizes[2] = vecDimNameSizes["lon"];
535  }
536  else if( varGridDims != nullptr )
537  {
538  // Obtain rectilinear attributes (dimension sizes)
539  NcDim* dimGridRank = varGridDims->get_dim( 0 );
540  if( dimGridRank == NULL )
541  {
542  _EXCEPTIONT( "Variable \"grid_dims\" has no dimensions" );
543  }
544 
545  int gridrank = dimGridRank->size();
546  // now get the rank dimensions
547  int gridsizes[2];
548  varGridDims->get( &( gridsizes[0] ), dimGridRank->size() );
549 
550  const Range& elems = remapper->GetMeshEntities( moab::Remapper::SourceMesh );
551  bool isQuads = elems.all_of_type( moab::MBQUAD );
552  bool isTris = elems.all_of_type( moab::MBTRI );
553 
554  if( gridrank == 1 )
555  {
556  vecDimSizes[0] = ( isQuads ? static_cast< int >( moab::TempestRemapper::CS )
557  : ( isTris ? static_cast< int >( moab::TempestRemapper::ICO )
558  : static_cast< int >( moab::TempestRemapper::ICOD ) ) );
559  vecDimSizes[1] = elems.size();
560  vecDimSizes[2] = remapper->GetMeshVertices( moab::Remapper::SourceMesh ).size();
561  }
562  else
563  {
564  vecDimSizes[0] = static_cast< int >( moab::TempestRemapper::RLL );
565  vecDimSizes[1] = gridsizes[0];
566  vecDimSizes[2] = gridsizes[1];
567  }
568  }
569  else
570  {
571  const Range& elems = remapper->GetMeshEntities( moab::Remapper::SourceMesh );
572  bool isQuads = elems.all_of_type( moab::MBQUAD );
573  bool isTris = elems.all_of_type( moab::MBTRI );
574  // vecDimSizes[0] = static_cast< int >( remapper->GetMeshType( moab::Remapper::SourceMesh );
575  vecDimSizes[0] = ( isQuads ? static_cast< int >( moab::TempestRemapper::CS )
576  : ( isTris ? static_cast< int >( moab::TempestRemapper::ICO )
577  : static_cast< int >( moab::TempestRemapper::ICOD ) ) );
578  vecDimSizes[1] = elems.size();
579  vecDimSizes[2] = remapper->GetMeshVertices( moab::Remapper::SourceMesh ).size();
580  }
581  moab::EntityHandle mSet = 0;
582  // mSet = remapper->GetMeshSet( moab::Remapper::SourceMesh );
583  result = gMB->tag_set_data( rectilinearTag, &mSet, 1, vecDimSizes.data() );
584  MB_CHK_ERR( result );
585 
586  switch( vecDimSizes[0] )
587  {
588  case 0:
589  printf( "Cubed-Sphere mesh: %d (elements), %d (vertices)\n", vecDimSizes[1], vecDimSizes[2] );
590  break;
591  case 1:
592  printf( "Rectilinear RLL mesh: (lon) %d X (lat) %d\n", vecDimSizes[2], vecDimSizes[1] );
593  break;
594  case 2:
595  printf( "Icosahedral (triangular) mesh: %d (elements), %d (vertices)\n", vecDimSizes[1],
596  vecDimSizes[2] );
597  break;
598  case 3:
599  default:
600  printf( "Polygonal mesh: %d (elements), %d (vertices)\n", vecDimSizes[1], vecDimSizes[2] );
601  break;
602  }
603 
604  ncInput.close();
605 
606  const size_t nOverlapFaces = tempestMesh->faces.size();
607  if( tempestMesh->vecSourceFaceIx.size() == nOverlapFaces &&
608  tempestMesh->vecSourceFaceIx.size() == nOverlapFaces )
609  {
610  int defaultInt = -1;
611  use_overlap_context = true;
612  // Check if our MOAB mesh has RED and BLUE tags; this would indicate we are
613  // converting an overlap grid
614  result = gMB->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag,
615  MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );
616  MB_CHK_SET_ERR( result, "can't create target parent tag" );
617 
618  result = gMB->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag,
619  MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );
620  MB_CHK_SET_ERR( result, "can't create source parent tag" );
621 
622  const Range& faces = remapper->GetMeshEntities( moab::Remapper::SourceMesh );
623 
624  std::vector< int > gids( faces.size() ), srcpar( faces.size() ), tgtpar( faces.size() );
625  result = gMB->tag_get_data( id_tag, faces, &gids[0] );
626  MB_CHK_ERR( result );
627 
628  for( unsigned ii = 0; ii < faces.size(); ++ii )
629  {
630  srcpar[ii] = tempestMesh->vecSourceFaceIx[gids[ii] - 1];
631  tgtpar[ii] = tempestMesh->vecTargetFaceIx[gids[ii] - 1];
632  }
633 
634  result = gMB->tag_set_data( srcParentTag, faces, &srcpar[0] );
635  MB_CHK_ERR( result );
636  result = gMB->tag_set_data( tgtParentTag, faces, &tgtpar[0] );
637  MB_CHK_ERR( result );
638 
639  srcpar.clear();
640  tgtpar.clear();
641  gids.clear();
642  }
643  }
644  else if( tempestout )
645  {
646  moab::EntityHandle& srcmesh = remapper->GetMeshSet( moab::Remapper::SourceMesh );
647  moab::EntityHandle& ovmesh = remapper->GetMeshSet( moab::Remapper::OverlapMesh );
648 
649  // load the mesh in MOAB format
650  std::vector< int > metadata( 2 );
651  result = remapper->LoadNativeMesh( *j, srcmesh, metadata );
652  MB_CHK_ERR( result );
653 
654  if( unitscaling )
655  {
656  result = moab::IntxUtils::ScaleToRadius( gMB, srcmesh, 1.0 );
657  MB_CHK_ERR( result );
658  }
659 
660  // Check if our MOAB mesh has RED and BLUE tags; this would indicate we are converting
661  // an overlap grid
662  ErrorCode rval1 = gMB->tag_get_handle( "SourceParent", srcParentTag );
663  ErrorCode rval2 = gMB->tag_get_handle( "TargetParent", tgtParentTag );
664  if( rval1 == MB_SUCCESS && rval2 == MB_SUCCESS )
665  {
666  use_overlap_context = true;
667  ovmesh = srcmesh;
668 
669  Tag countTag;
670  result = gMB->tag_get_handle( "Counting", countTag );
671  MB_CHK_ERR( result );
672 
673  // Load the meshes and validate
674  Tag order;
675  ReorderTool reorder_tool( &core );
676  result = reorder_tool.handle_order_from_int_tag( srcParentTag, -1, order );
677  MB_CHK_ERR( result );
678  result = reorder_tool.reorder_entities( order );
679  MB_CHK_ERR( result );
680  result = gMB->tag_delete( order );
681  MB_CHK_ERR( result );
682  result = remapper->ConvertMeshToTempest( moab::Remapper::OverlapMesh );
683  MB_CHK_ERR( result );
684  }
685  else
686  {
687  if( metadata[0] == static_cast< int >( moab::TempestRemapper::RLL ) )
688  {
689  assert( metadata.size() );
690  std::cout << "Converting a RLL mesh with rectilinear dimension: " << metadata[0] << " X "
691  << metadata[1] << std::endl;
692  }
693 
694  // Convert the mesh and validate
695  result = remapper->ConvertMeshToTempest( moab::Remapper::SourceMesh );
696  MB_CHK_ERR( result );
697  }
698  }
699  }
700  else
701  result = gMB->load_file( j->c_str(), 0, read_options.c_str() );
702 #else
703  result = gMB->load_file( j->c_str(), 0, read_options.c_str() );
704 #endif
705  if( MB_SUCCESS != result )
706  {
707  std::cerr << "Failed to load \"" << *j << "\"." << std::endl;
708  std::cerr << "Error code: " << gMB->get_error_string( result ) << " (" << result << ")" << std::endl;
709  std::string message;
710  if( MB_SUCCESS == gMB->get_last_error( message ) && !message.empty() )
711  std::cerr << "Error message: " << message << std::endl;
712 #ifdef MOAB_HAVE_MPI
713  MPI_Finalize();
714 #endif
715  return READ_ERROR;
716  }
717  if( !proc_id ) std::cerr << "Read \"" << *j << "\"" << std::endl;
718  if( print_times && !proc_id ) write_times( std::cout );
719  }
720 
721  // Determine if the user has specified any geometry sets to write
722  bool have_geom = false;
723  for( dim = 0; dim <= 3; ++dim )
724  {
725  if( !geom[dim].empty() ) have_geom = true;
726  if( verbose ) print_id_list( geom_names[dim], std::cout, geom[dim] );
727  }
728 
729  // True if the user has specified any sets to write
730  bool have_sets = have_geom;
731 
732  // Get geometry tags
733  Tag dim_tag;
734  if( have_geom )
735  {
736  if( id_tag == 0 )
737  {
738  std::cerr << "No ID tag defined." << std::endl;
739  have_geom = false;
740  }
741  result = gMB->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dim_tag );
742  if( MB_SUCCESS != result )
743  {
744  std::cerr << "No geometry tag defined." << std::endl;
745  have_geom = false;
746  }
747  }
748 
749  // Get geometry sets
750  if( have_geom )
751  {
752  int id_val;
753  Tag tags[] = { id_tag, dim_tag };
754  const void* vals[] = { &id_val, &dim };
755  for( dim = 0; dim <= 3; ++dim )
756  {
757  int init_count = set_list.size();
758  for( std::set< int >::iterator iter = geom[dim].begin(); iter != geom[dim].end(); ++iter )
759  {
760  id_val = *iter;
761  range.clear();
762  result = gMB->get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, range );
763  if( MB_SUCCESS != result || range.empty() )
764  {
765  range.clear();
766  std::cerr << geom_names[dim] << " " << id_val << " not found.\n";
767  }
768  std::copy( range.begin(), range.end(), std::back_inserter( set_list ) );
769  }
770 
771  if( verbose )
772  std::cout << "Found " << ( set_list.size() - init_count ) << ' ' << geom_names[dim] << " sets"
773  << std::endl;
774  }
775  }
776 
777  // Get mesh groupings
778  for( i = 0; i < 4; ++i )
779  {
780  if( verbose ) print_id_list( mesh_tag_names[i], std::cout, mesh[i] );
781 
782  if( mesh[i].empty() ) continue;
783  have_sets = true;
784 
785  // Get tag
786  Tag tag;
787  result = gMB->tag_get_handle( mesh_tag_names[i], 1, MB_TYPE_INTEGER, tag );
788  if( MB_SUCCESS != result )
789  {
790  std::cerr << "Tag not found: " << mesh_tag_names[i] << std::endl;
791  continue;
792  }
793 
794  // get entity sets
795  int init_count = set_list.size();
796  for( std::set< int >::iterator iter = mesh[i].begin(); iter != mesh[i].end(); ++iter )
797  {
798  range.clear();
799  const void* vals[] = { &*iter };
800  result = gMB->get_entities_by_type_and_tag( 0, MBENTITYSET, &tag, vals, 1, range );
801  if( MB_SUCCESS != result || range.empty() )
802  {
803  range.clear();
804  std::cerr << mesh_tag_names[i] << " " << *iter << " not found.\n";
805  }
806  std::copy( range.begin(), range.end(), std::back_inserter( set_list ) );
807  }
808 
809  if( verbose )
810  std::cout << "Found " << ( set_list.size() - init_count ) << ' ' << mesh_tag_names[i] << " sets"
811  << std::endl;
812  }
813 
814  // Check if output is limited to certain dimensions of elements
815  bool bydim = false;
816  for( dim = 1; dim < 4; ++dim )
817  if( dims[dim] ) bydim = true;
818 
819  // Check conflicting input
820  if( bydim )
821  {
822  if( generate[1] && !dims[1] )
823  {
824  std::cerr << "Warning: Request to generate 1D internal entities but not export them." << std::endl;
825  generate[1] = false;
826  }
827  if( generate[2] && !dims[2] )
828  {
829  std::cerr << "Warning: Request to generate 2D internal entities but not export them." << std::endl;
830  generate[2] = false;
831  }
832  }
833 
834  // Generate any internal entities
835  if( generate[1] || generate[2] )
836  {
837  EntityHandle all_mesh = 0;
838  const EntityHandle* sets = &all_mesh;
839  int num_sets = 1;
840  if( have_sets )
841  {
842  num_sets = set_list.size();
843  sets = &set_list[0];
844  }
845  for( i = 0; i < num_sets; ++i )
846  {
847  Range dim3, dim2, adj;
848  gMB->get_entities_by_dimension( sets[i], 3, dim3, true );
849  if( generate[1] )
850  {
851  gMB->get_entities_by_dimension( sets[i], 2, dim2, true );
852  gMB->get_adjacencies( dim3, 1, true, adj, Interface::UNION );
853  gMB->get_adjacencies( dim2, 1, true, adj, Interface::UNION );
854  }
855  if( generate[2] )
856  {
857  gMB->get_adjacencies( dim3, 2, true, adj, Interface::UNION );
858  }
859  if( sets[i] ) gMB->add_entities( sets[i], adj );
860  }
861  }
862 
863  // Delete any entities not of the dimensions to be exported
864  if( bydim )
865  {
866  // Get list of dead elements
867  Range dead_entities, tmp_range;
868  for( dim = 1; dim <= 3; ++dim )
869  {
870  if( dims[dim] ) continue;
871  gMB->get_entities_by_dimension( 0, dim, tmp_range );
872  dead_entities.merge( tmp_range );
873  }
874  // Remove dead entities from all sets, and add all
875  // empty sets to list of dead entities.
876  Range empty_sets;
877  remove_entities_from_sets( gMB, dead_entities, empty_sets );
878  while( !empty_sets.empty() )
879  {
880  if( !set_list.empty() ) remove_from_vector( set_list, empty_sets );
881  dead_entities.merge( empty_sets );
882  tmp_range.clear();
883  remove_entities_from_sets( gMB, empty_sets, tmp_range );
884  empty_sets = subtract( tmp_range, dead_entities );
885  }
886  // Destroy dead entities
887  gMB->delete_entities( dead_entities );
888  }
889 
890  // If user specified sets to write, but none were found, exit.
891  if( have_sets && set_list.empty() )
892  {
893  std::cerr << "Nothing to write." << std::endl;
894 #ifdef MOAB_HAVE_MPI
895  MPI_Finalize();
896 #endif
897  return ENT_NOT_FOUND;
898  }
899 
900  // interpret the mpas partition file created by gpmetis
901  if( !metis_partition_file.empty() )
902  {
903  int err = process_partition_file( gMB, metis_partition_file );
904  if( err )
905  {
906  std::cerr << "Failed to process partition file \"" << metis_partition_file << "\"." << std::endl;
907 #ifdef MOAB_HAVE_MPI
908  MPI_Finalize();
909 #endif
910  return WRITE_ERROR;
911  }
912  }
913  if( verbose )
914  {
915  if( have_sets )
916  std::cout << "Found " << set_list.size() << " specified sets to write (total)." << std::endl;
917  else
918  std::cout << "No sets specifed. Writing entire mesh." << std::endl;
919  }
920 
921  // Write the output file
922  reset_times();
923 #ifdef MOAB_HAVE_TEMPESTREMAP
924  if( remapper )
925  {
926  Range faces;
927  Mesh* tempestMesh =
928  remapper->GetMesh( ( use_overlap_context ? moab::Remapper::OverlapMesh : moab::Remapper::SourceMesh ) );
929  moab::EntityHandle& srcmesh =
930  remapper->GetMeshSet( ( use_overlap_context ? moab::Remapper::OverlapMesh : moab::Remapper::SourceMesh ) );
931  result = gMB->get_entities_by_dimension( srcmesh, 2, faces );
932  MB_CHK_ERR( result );
933  int ntot_elements = 0, nelements = faces.size();
934 #ifdef MOAB_HAVE_MPI
935  int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, pcomm->comm() );
936  if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
937 #else
938  ntot_elements = nelements;
939 #endif
940 
941  Tag gidTag = gMB->globalId_tag();
942  std::vector< int > gids( faces.size() );
943  result = gMB->tag_get_data( gidTag, faces, &gids[0] );
944  MB_CHK_ERR( result );
945 
946  if( faces.size() > 1 && gids[0] == gids[1] && !use_overlap_context )
947  {
948 #ifdef MOAB_HAVE_MPI
949  result = pcomm->assign_global_ids( srcmesh, 2, 1, false );
950  MB_CHK_ERR( result );
951 #else
952  result = remapper->assign_vertex_element_IDs( gidTag, srcmesh, 2, 1 );
953  MB_CHK_ERR( result );
954  result = remapper->assign_vertex_element_IDs( gidTag, srcmesh, 0, 1 );
955  MB_CHK_ERR( result );
956 #endif
957  }
958 
959  // VSM: If user requested explicitly for some metadata, we need to generate the DoF ID tag
960  // and set the appropriate numbering based on specified discretization order
961  // Useful only for SE meshes with GLL DoFs
962  if( spectral_order > 1 && globalid_tag_name.size() > 1 )
963  {
964  result = remapper->GenerateMeshMetadata( *tempestMesh, ntot_elements, faces, NULL, globalid_tag_name,
965  spectral_order );
966  MB_CHK_ERR( result );
967  }
968 
969  if( tempestout )
970  {
971  // Check if our MOAB mesh has RED and BLUE tags; this would indicate we are converting an
972  // overlap grid
973  if( use_overlap_context && false )
974  {
975  const int nOverlapFaces = faces.size();
976  // Overlap mesh: resize the source and target connection arrays
977  tempestMesh->vecSourceFaceIx.resize( nOverlapFaces ); // 0-based indices corresponding to source mesh
978  tempestMesh->vecTargetFaceIx.resize( nOverlapFaces ); // 0-based indices corresponding to target mesh
979  result = gMB->tag_get_data( srcParentTag, faces, &tempestMesh->vecSourceFaceIx[0] );
980  MB_CHK_ERR( result );
981  result = gMB->tag_get_data( tgtParentTag, faces, &tempestMesh->vecTargetFaceIx[0] );
982  MB_CHK_ERR( result );
983  }
984  // Write out the mesh using TempestRemap
985  tempestMesh->Write( out, NcFile::Netcdf4 );
986  file_written = true;
987  }
988  delete remapper; // cleanup
989  }
990 #endif
991 
992  if( !file_written )
993  {
994  if( have_sets )
995  result = gMB->write_file( out.c_str(), format, write_options.c_str(), &set_list[0], set_list.size() );
996  else
997  result = gMB->write_file( out.c_str(), format, write_options.c_str() );
998  if( MB_SUCCESS != result )
999  {
1000  std::cerr << "Failed to write \"" << out << "\"." << std::endl;
1001  std::cerr << "Error code: " << gMB->get_error_string( result ) << " (" << result << ")" << std::endl;
1002  std::string message;
1003  if( MB_SUCCESS == gMB->get_last_error( message ) && !message.empty() )
1004  std::cerr << "Error message: " << message << std::endl;
1005 #ifdef MOAB_HAVE_MPI
1006  MPI_Finalize();
1007 #endif
1008  return WRITE_ERROR;
1009  }
1010  }
1011 
1012  if( !proc_id ) std::cerr << "Wrote \"" << out << "\"" << std::endl;
1013  if( print_times && !proc_id ) write_times( std::cout );
1014 
1015 #ifdef MOAB_HAVE_MPI
1016  MPI_Finalize();
1017 #endif
1018  return 0;
1019 }
1020 
1021 bool parse_id_list( const char* string, std::set< int >& results )
1022 {
1023  bool okay = true;
1024  char* mystr = strdup( string );
1025  for( const char* ptr = strtok( mystr, "," ); ptr; ptr = strtok( 0, "," ) )
1026  {
1027  char* endptr;
1028  long val = strtol( ptr, &endptr, 0 );
1029  if( endptr == ptr || val <= 0 )
1030  {
1031  std::cerr << "Not a valid id: " << ptr << std::endl;
1032  okay = false;
1033  break;
1034  }
1035 
1036  long val2 = val;
1037  if( *endptr == '-' )
1038  {
1039  const char* sptr = endptr + 1;
1040  val2 = strtol( sptr, &endptr, 0 );
1041  if( endptr == sptr || val2 <= 0 )
1042  {
1043  std::cerr << "Not a valid id: " << sptr << std::endl;
1044  okay = false;
1045  break;
1046  }
1047  if( val2 < val )
1048  {
1049  std::cerr << "Invalid id range: " << ptr << std::endl;
1050  okay = false;
1051  break;
1052  }
1053  }
1054 
1055  if( *endptr )
1056  {
1057  std::cerr << "Unexpected character: " << *endptr << std::endl;
1058  okay = false;
1059  break;
1060  }
1061 
1062  for( ; val <= val2; ++val )
1063  if( !results.insert( (int)val ).second ) std::cerr << "Warning: duplicate Id: " << val << std::endl;
1064  }
1065 
1066  free( mystr );
1067  return okay;
1068 }
1069 
1070 void print_id_list( const char* head, std::ostream& stream, const std::set< int >& list )
1071 {
1072  stream << head << ": ";
1073 
1074  if( list.empty() )
1075  {
1076  stream << "(none)" << std::endl;
1077  return;
1078  }
1079 
1080  int start, prev;
1081  std::set< int >::const_iterator iter = list.begin();
1082  start = prev = *( iter++ );
1083  for( ;; )
1084  {
1085  if( iter == list.end() || *iter != 1 + prev )
1086  {
1087  stream << start;
1088  if( prev != start ) stream << '-' << prev;
1089  if( iter == list.end() ) break;
1090  stream << ", ";
1091  start = *iter;
1092  }
1093  prev = *( iter++ );
1094  }
1095 
1096  stream << std::endl;
1097 }
1098 
1099 static void print_time( int clk_per_sec, const char* prefix, clock_t ticks, std::ostream& stream )
1100 {
1101  ticks *= clk_per_sec / 100;
1102  clock_t centi = ticks % 100;
1103  clock_t seconds = ticks / 100;
1104  stream << prefix;
1105  if( seconds < 120 )
1106  {
1107  stream << ( ticks / 100 ) << "." << centi << "s" << std::endl;
1108  }
1109  else
1110  {
1111  clock_t minutes = ( seconds / 60 ) % 60;
1112  clock_t hours = ( seconds / 3600 );
1113  seconds %= 60;
1114  if( hours ) stream << hours << "h";
1115  if( minutes ) stream << minutes << "m";
1116  if( seconds || centi ) stream << seconds << "." << centi << "s";
1117  stream << " (" << ( ticks / 100 ) << "." << centi << "s)" << std::endl;
1118  }
1119 }
1120 
1121 clock_t usr_time, sys_time, abs_time;
1122 
1123 #ifdef WIN32
1124 
1125 void reset_times()
1126 {
1127  abs_time = clock();
1128 }
1129 
1130 void write_times( std::ostream& stream )
1131 {
1132  clock_t abs_tm = clock();
1133  print_time( CLOCKS_PER_SEC, " ", abs_tm - abs_time, stream );
1134  abs_time = abs_tm;
1135 }
1136 
1137 #else
1138 
1139 void reset_times()
1140 {
1141  tms timebuf;
1142  abs_time = times( &timebuf );
1143  usr_time = timebuf.tms_utime;
1144  sys_time = timebuf.tms_stime;
1145 }
1146 
1147 void write_times( std::ostream& stream )
1148 {
1149  clock_t usr_tm, sys_tm, abs_tm;
1150  tms timebuf;
1151  abs_tm = times( &timebuf );
1152  usr_tm = timebuf.tms_utime;
1153  sys_tm = timebuf.tms_stime;
1154  print_time( sysconf( _SC_CLK_TCK ), " real: ", abs_tm - abs_time, stream );
1155  print_time( sysconf( _SC_CLK_TCK ), " user: ", usr_tm - usr_time, stream );
1156  print_time( sysconf( _SC_CLK_TCK ), " system: ", sys_tm - sys_time, stream );
1157  abs_time = abs_tm;
1158  usr_time = usr_tm;
1159  sys_time = sys_tm;
1160 }
1161 
1162 #endif
1163 
1164 bool make_opts_string( std::vector< std::string > options, std::string& opts )
1165 {
1166  opts.clear();
1167  if( options.empty() ) return true;
1168 
1169  // choose a separator character
1170  std::vector< std::string >::const_iterator i;
1171  char separator = '\0';
1172  const char* alt_separators = ";+,:\t\n";
1173  for( const char* sep_ptr = alt_separators; *sep_ptr; ++sep_ptr )
1174  {
1175  bool seen = false;
1176  for( i = options.begin(); i != options.end(); ++i )
1177  if( i->find( *sep_ptr, 0 ) != std::string::npos )
1178  {
1179  seen = true;
1180  break;
1181  }
1182  if( !seen )
1183  {
1184  separator = *sep_ptr;
1185  break;
1186  }
1187  }
1188  if( !separator )
1189  {
1190  std::cerr << "Error: cannot find separator character for options string" << std::endl;
1191  return false;
1192  }
1193  if( separator != ';' )
1194  {
1195  opts = ";";
1196  opts += separator;
1197  }
1198 
1199  // concatenate options
1200  i = options.begin();
1201  opts += *i;
1202  for( ++i; i != options.end(); ++i )
1203  {
1204  opts += separator;
1205  opts += *i;
1206  }
1207 
1208  return true;
1209 }
1210 
1211 void list_formats( Interface* gMB )
1212 {
1213  const char iface_name[] = "ReaderWriterSet";
1214  ErrorCode err;
1215  ReaderWriterSet* set = 0;
1216  ReaderWriterSet::iterator i;
1217  std::ostream& str = std::cout;
1218 
1219  // get ReaderWriterSet
1220  err = gMB->query_interface( set );
1221  if( err != MB_SUCCESS || !set )
1222  {
1223  std::cerr << "Internal error: Interface \"" << iface_name << "\" not available.\n";
1224  exit( OTHER_ERROR );
1225  }
1226 
1227  // get field width for format description
1228  size_t w = 0;
1229  for( i = set->begin(); i != set->end(); ++i )
1230  if( i->description().length() > w ) w = i->description().length();
1231 
1232  // write table header
1233  str << "Format " << std::setw( w ) << std::left << "Description" << " Read Write File Name Suffixes\n"
1234  << "------ " << std::setw( w ) << std::setfill( '-' ) << "" << std::setfill( ' ' )
1235  << " ---- ----- ------------------\n";
1236 
1237  // write table data
1238  for( i = set->begin(); i != set->end(); ++i )
1239  {
1240  std::vector< std::string > ext;
1241  i->get_extensions( ext );
1242  str << std::setw( 6 ) << i->name() << " " << std::setw( w ) << std::left << i->description() << " "
1243  << ( i->have_reader() ? " yes" : " no" ) << " " << ( i->have_writer() ? " yes" : " no" ) << " ";
1244  for( std::vector< std::string >::iterator j = ext.begin(); j != ext.end(); ++j )
1245  str << " " << *j;
1246  str << std::endl;
1247  }
1248  str << std::endl;
1249 
1250  gMB->release_interface( set );
1251  exit( 0 );
1252 }
1253 
1254 void remove_entities_from_sets( Interface* gMB, Range& dead_entities, Range& empty_sets )
1255 {
1256  empty_sets.clear();
1257  Range sets;
1258  gMB->get_entities_by_type( 0, MBENTITYSET, sets );
1259  for( Range::iterator i = sets.begin(); i != sets.end(); ++i )
1260  {
1261  Range set_contents;
1262  gMB->get_entities_by_handle( *i, set_contents, false );
1263  set_contents = intersect( set_contents, dead_entities );
1264  gMB->remove_entities( *i, set_contents );
1265  set_contents.clear();
1266  gMB->get_entities_by_handle( *i, set_contents, false );
1267  if( set_contents.empty() ) empty_sets.insert( *i );
1268  }
1269 }
1270 
1271 void remove_from_vector( std::vector< EntityHandle >& vect, const Range& ents_to_remove )
1272 {
1273  Range::const_iterator i;
1274  std::vector< EntityHandle >::iterator j;
1275  for( i = ents_to_remove.begin(); i != ents_to_remove.end(); ++i )
1276  {
1277  j = std::find( vect.begin(), vect.end(), *i );
1278  if( j != vect.end() ) vect.erase( j );
1279  }
1280 }
1281 
1282 std::string percent_subst( const std::string& s, int val )
1283 {
1284  if( s.empty() ) return s;
1285 
1286  size_t j = s.find( '%' );
1287  if( j == std::string::npos ) return s;
1288 
1289  std::ostringstream st;
1290  st << s.substr( 0, j );
1291  st << val;
1292 
1293  size_t i;
1294  while( ( i = s.find( '%', j + 1 ) ) != std::string::npos )
1295  {
1296  st << s.substr( j, i - j );
1297  st << val;
1298  j = i;
1299  }
1300  st << s.substr( j + 1 );
1301  return st.str();
1302 }
1303 
1304 int process_partition_file( Interface* mb, std::string& metis_partition_file )
1305 {
1306  // how many faces in the file ? how do we make sure it is an mpas file?
1307  // mpas atmosphere files can be downloaded from here
1308  // https://mpas-dev.github.io/atmosphere/atmosphere_meshes.html
1309  Range faces;
1310  MB_CHK_ERR( mb->get_entities_by_dimension( 0, 2, faces ) );
1311  std::cout << " MPAS model has " << faces.size() << " polygons\n";
1312 
1313  // read the partition file
1314  std::ifstream partfile;
1315  partfile.open( metis_partition_file.c_str() );
1316  std::string line;
1317  std::vector< int > parts;
1318  parts.resize( faces.size(), -1 );
1319  int i = 0;
1320  if( partfile.is_open() )
1321  {
1322  while( getline( partfile, line ) )
1323  {
1324  // cout << line << '\n';
1325  parts[i++] = atoi( line.c_str() );
1326  if( i > (int)faces.size() )
1327  {
1328  std::cerr << " too many lines in partition file \n. bail out \n";
1329  return 1;
1330  }
1331  }
1332  partfile.close();
1333  }
1334  std::vector< int >::iterator pmax = max_element( parts.begin(), parts.end() );
1335  std::vector< int >::iterator pmin = min_element( parts.begin(), parts.end() );
1336  if( *pmin <= -1 )
1337  {
1338  std::cerr << " partition file is incomplete, *pmin is -1 !! \n";
1339  return 1;
1340  }
1341  std::cout << " partitions range: " << *pmin << " " << *pmax << "\n";
1342  Tag part_set_tag;
1343  int dum_id = -1;
1344  MB_CHK_ERR( mb->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_set_tag,
1345  MB_TAG_SPARSE | MB_TAG_CREAT, &dum_id ) );
1346 
1347  // get any sets already with this tag, and clear them
1348  // remove the parallel partition sets if they exist
1349  Range tagged_sets;
1350  MB_CHK_ERR(
1351  mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &part_set_tag, NULL, 1, tagged_sets, Interface::UNION ) );
1352  if( !tagged_sets.empty() )
1353  {
1354  MB_CHK_ERR( mb->clear_meshset( tagged_sets ) );
1355  MB_CHK_ERR( mb->tag_delete_data( part_set_tag, tagged_sets ) );
1356  }
1357  Tag gid;
1358  MB_CHK_ERR( mb->tag_get_handle( "GLOBAL_ID", gid ) );
1359  int num_sets = *pmax + 1;
1360  if( *pmin != 0 )
1361  {
1362  std::cout << " problem reading parts; min is not 0 \n";
1363  return 1;
1364  }
1365  for( i = 0; i < num_sets; i++ )
1366  {
1367  EntityHandle new_set;
1368  MB_CHK_ERR( mb->create_meshset( MESHSET_SET, new_set ) );
1369  tagged_sets.insert( new_set );
1370  }
1371  int* dum_ids = new int[num_sets];
1372  for( i = 0; i < num_sets; i++ )
1373  dum_ids[i] = i;
1374 
1375  MB_CHK_ERR( mb->tag_set_data( part_set_tag, tagged_sets, dum_ids ) );
1376  delete[] dum_ids;
1377 
1378  std::vector< int > gids;
1379  int num_faces = (int)faces.size();
1380  gids.resize( num_faces );
1381  MB_CHK_ERR( mb->tag_get_data( gid, faces, &gids[0] ) );
1382 
1383  for( int j = 0; j < num_faces; j++ )
1384  {
1385  int eid = gids[j];
1386  EntityHandle eh = faces[j];
1387  int partition = parts[eid - 1];
1388  if( partition < 0 || partition >= num_sets )
1389  {
1390  std::cout << " wrong partition number \n";
1391  return 1;
1392  }
1393  MB_CHK_ERR( mb->add_entities( tagged_sets[partition], &eh, 1 ) );
1394  }
1395  return 0;
1396 }