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