Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
NCHelperDomain.cpp
Go to the documentation of this file.
1 #include "NCHelperDomain.hpp"
2 #include "moab/FileOptions.hpp"
3 #include "moab/ReadUtilIface.hpp"
5 #include "AEntityFactory.hpp"
6 #ifdef MOAB_HAVE_MPI
8 #endif
9 #include <cmath>
10 #include <sstream>
11 
12 namespace moab
13 {
14 
15 bool NCHelperDomain::can_read_file( ReadNC* readNC, int fileId )
16 {
17  std::vector< std::string >& dimNames = readNC->dimNames;
18 
19  // If dimension names "n" AND "ni" AND "nj" AND "nv" exist then it should be the Domain grid
20  if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "n" ) ) != dimNames.end() ) &&
21  ( std::find( dimNames.begin(), dimNames.end(), std::string( "ni" ) ) != dimNames.end() ) &&
22  ( std::find( dimNames.begin(), dimNames.end(), std::string( "nj" ) ) != dimNames.end() ) &&
23  ( std::find( dimNames.begin(), dimNames.end(), std::string( "nv" ) ) != dimNames.end() ) )
24  {
25  // Make sure it is CAM grid
26  std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "source" );
27  if( attIt == readNC->globalAtts.end() ) return false;
28  unsigned int sz = attIt->second.attLen;
29  std::string att_data;
30  att_data.resize( sz + 1 );
31  att_data[sz] = '\000';
32  int success =
33  NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
34  if( success ) return false;
35  /*if (att_data.find("CAM") == std::string::npos)
36  return false;*/
37 
38  return true;
39  }
40 
41  return false;
42 }
43 
45 {
46  Interface*& mbImpl = _readNC->mbImpl;
47  std::vector< std::string >& dimNames = _readNC->dimNames;
48  std::vector< int >& dimLens = _readNC->dimLens;
49  std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
50  DebugOutput& dbgOut = _readNC->dbgOut;
51 #ifdef MOAB_HAVE_MPI
52  bool& isParallel = _readNC->isParallel;
53 #endif
54  int& partMethod = _readNC->partMethod;
55  ScdParData& parData = _readNC->parData;
56 
57  ErrorCode rval;
58 
59  // Look for names of i/j dimensions
60  // First i
61  std::vector< std::string >::iterator vit;
62  unsigned int idx;
63  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ni" ) ) != dimNames.end() )
64  idx = vit - dimNames.begin();
65  else
66  {
67  MB_SET_ERR( MB_FAILURE, "Couldn't find 'ni' variable" );
68  }
69  iDim = idx;
70  gDims[0] = 0;
71  gDims[3] = dimLens[idx];
72 
73  // Then j
74  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "nj" ) ) != dimNames.end() )
75  idx = vit - dimNames.begin();
76  else
77  {
78  MB_SET_ERR( MB_FAILURE, "Couldn't find 'nj' variable" );
79  }
80  jDim = idx;
81  gDims[1] = 0;
82  gDims[4] = dimLens[idx]; // Add 2 for the pole points ? not needed
83 
84  // do not use gcdims ? or use only gcdims?
85 
86  // Try a truly 2D mesh
87  gDims[2] = -1;
88  gDims[5] = -1;
89 
90  // Get number of vertices per cell
91  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "nv" ) ) != dimNames.end() )
92  idx = vit - dimNames.begin();
93  else
94  {
95  MB_SET_ERR( MB_FAILURE, "Couldn't find 'nv' dimension" );
96  }
97  nvDim = idx;
98  nv = dimLens[idx];
99 
100  // Parse options to get subset
101  int rank = 0, procs = 1;
102 #ifdef MOAB_HAVE_MPI
103  if( isParallel )
104  {
105  ParallelComm*& myPcomm = _readNC->myPcomm;
106  rank = myPcomm->proc_config().proc_rank();
107  procs = myPcomm->proc_config().proc_size();
108  }
109 #endif
110  if( procs > 1 )
111  {
112  for( int i = 0; i < 6; i++ )
113  parData.gDims[i] = gDims[i];
114  parData.partMethod = partMethod;
115  int pdims[3];
116 
118  rval = ScdInterface::compute_partition( procs, rank, parData, lDims, locallyPeriodic, pdims );MB_CHK_ERR( rval );
119  for( int i = 0; i < 3; i++ )
120  parData.pDims[i] = pdims[i];
121 
122  dbgOut.tprintf( 1, "Partition: %dx%d (out of %dx%d)\n", lDims[3] - lDims[0], lDims[4] - lDims[1],
123  gDims[3] - gDims[0], gDims[4] - gDims[1] );
124  if( 0 == rank )
125  dbgOut.tprintf( 1, "Contiguous chunks of size %d bytes.\n",
126  8 * ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ) );
127  }
128  else
129  {
130  for( int i = 0; i < 6; i++ )
131  lDims[i] = gDims[i];
133  }
134 
135  // Now get actual coordinate values for vertices and cell centers
136  lCDims[0] = lDims[0];
137 
138  lCDims[3] = lDims[3];
139 
140  // For FV models, will always be non-periodic in j
141  lCDims[1] = lDims[1];
142  lCDims[4] = lDims[4];
143 
144 #if 0
145  // Resize vectors to store values later
146  if (-1 != lDims[0])
147  ilVals.resize(lDims[3] - lDims[0] + 1);
148  if (-1 != lCDims[0])
149  ilCVals.resize(lCDims[3] - lCDims[0] + 1);
150  if (-1 != lDims[1])
151  jlVals.resize(lDims[4] - lDims[1] + 1);
152  if (-1 != lCDims[1])
153  jlCVals.resize(lCDims[4] - lCDims[1] + 1);
154  if (nTimeSteps > 0)
155  tVals.resize(nTimeSteps);
156 #endif
157 
158  dbgOut.tprintf( 1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4] );
159  dbgOut.tprintf( 1, "%d elements, %d vertices\n", ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ),
160  ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ) * nv );
161 
162  // For each variable, determine the entity location type and number of levels
163  std::map< std::string, ReadNC::VarData >::iterator mit;
164  for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
165  {
166  ReadNC::VarData& vd = ( *mit ).second;
167 
168  // Default entLoc is ENTLOCSET
169  if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
170  {
171  if( ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) &&
172  ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) )
174  else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jDim ) != vd.varDims.end() ) &&
175  ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) )
177  else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) &&
178  ( std::find( vd.varDims.begin(), vd.varDims.end(), iDim ) != vd.varDims.end() ) )
180  }
181 
182  // Default numLev is 0
183  if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels;
184  }
185 
186  std::vector< std::string > ijdimNames( 2 );
187  ijdimNames[0] = "__ni";
188  ijdimNames[1] = "__nj";
189  std::string tag_name;
190  Tag tagh;
191 
192  // __<dim_name>_LOC_MINMAX (for slon, slat, lon and lat)
193  for( unsigned int i = 0; i != ijdimNames.size(); i++ )
194  {
195  std::vector< int > val( 2, 0 );
196  if( ijdimNames[i] == "__ni" )
197  {
198  val[0] = lDims[0];
199  val[1] = lDims[3];
200  }
201  else if( ijdimNames[i] == "__nj" )
202  {
203  val[0] = lDims[1];
204  val[1] = lDims[4];
205  }
206 
207  std::stringstream ss_tag_name;
208  ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
209  tag_name = ss_tag_name.str();
210  rval = mbImpl->tag_get_handle( tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
211  rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
212  if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
213  }
214 
215  // __<dim_name>_LOC_VALS (for slon, slat, lon and lat)
216  // Assume all have the same data type as lon (expected type is float or double)
217  switch( varInfo["xc"].varDataType )
218  {
219  case NC_FLOAT:
220  case NC_DOUBLE:
221  break;
222  default:
223  MB_SET_ERR( MB_FAILURE, "Unexpected variable data type for 'lon'" );
224  }
225 
226  // do not need conventional tags
227  Tag convTagsCreated = 0;
228  int def_val = 0;
229  rval = mbImpl->tag_get_handle( "__CONV_TAGS_CREATED", 1, MB_TYPE_INTEGER, convTagsCreated,
230  MB_TAG_SPARSE | MB_TAG_CREAT, &def_val );MB_CHK_SET_ERR( rval, "Trouble getting _CONV_TAGS_CREATED tag" );
231  int create_conv_tags_flag = 1;
232  rval = mbImpl->tag_set_data( convTagsCreated, &_fileSet, 1, &create_conv_tags_flag );MB_CHK_SET_ERR( rval, "Trouble setting _CONV_TAGS_CREATED tag" );
233 
234  return MB_SUCCESS;
235 }
236 
238 {
239  Interface*& mbImpl = _readNC->mbImpl;
240  // std::string& fileName = _readNC->fileName;
241  Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
242  // const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
243  DebugOutput& dbgOut = _readNC->dbgOut;
244  /*int& gatherSetRank = _readNC->gatherSetRank;
245  int& trivialPartitionShift = _readNC->trivialPartitionShift;*/
246  /*
247 
248  int rank = 0;
249  int procs = 1;
250  #ifdef MOAB_HAVE_MPI
251  bool& isParallel = _readNC->isParallel;
252  if (isParallel) {
253  ParallelComm*& myPcomm = _readNC->myPcomm;
254  rank = myPcomm->proc_config().proc_rank();
255  procs = myPcomm->proc_config().proc_size();
256  }
257  #endif
258  */
259 
260  ErrorCode rval;
261  int success = 0;
262 
263  /*
264  bool create_gathers = false;
265  if (rank == gatherSetRank)
266  create_gathers = true;
267 
268  // Shift rank to obtain a rotated trivial partition
269  int shifted_rank = rank;
270  if (procs >= 2 && trivialPartitionShift > 0)
271  shifted_rank = (rank + trivialPartitionShift) % procs;*/
272 
273  // how many will have mask 0 or 1
274  // how many will have a fraction ? we will not instantiate all elements; only those with mask 1
275  // ? also, not all vertices, only those that belong to mask 1 elements ? we will not care about
276  // duplicate vertices; maybe another time ? we will start reading masks, vertices
277  int local_elems = ( lDims[4] - lDims[1] ) * ( lDims[3] - lDims[0] );
278  dbgOut.tprintf( 1, "local cells: %d \n", local_elems );
279 
280  // count how many will be with mask 1 here
281  // basically, read the mask variable on the local elements;
282  std::string maskstr( "mask" );
283  ReadNC::VarData& vmask = _readNC->varInfo[maskstr];
284 
285  // mask is (nj, ni)
286  vmask.readStarts.push_back( lDims[1] );
287  vmask.readStarts.push_back( lDims[0] );
288  vmask.readCounts.push_back( lDims[4] - lDims[1] );
289  vmask.readCounts.push_back( lDims[3] - lDims[0] );
290  std::vector< int > mask( local_elems );
291  success = NCFUNCAG( _vara_int )( _fileId, vmask.varId, &vmask.readStarts[0], &vmask.readCounts[0], &mask[0] );
292  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read int data for mask variable " );
293 
294  int nb_with_mask1 = 0;
295  for( int i = 0; i < local_elems; i++ )
296  if( 1 == mask[i] ) nb_with_mask1++;
297 
298  dbgOut.tprintf( 1, "local cells with mask 1: %d \n", nb_with_mask1 );
299  std::vector< NCDF_SIZE > startsv( 3 );
300  startsv[0] = vmask.readStarts[0];
301  startsv[1] = vmask.readStarts[1];
302  startsv[2] = 0;
303  std::vector< NCDF_SIZE > countsv( 3 );
304  countsv[0] = vmask.readCounts[0];
305  countsv[1] = vmask.readCounts[1];
306  countsv[2] = nv; // number of vertices per element
307 
308  // read xv and yv coords for vertices, and create elements;
309  std::string xvstr( "xv" );
310  ReadNC::VarData& var_xv = _readNC->varInfo[xvstr];
311  std::vector< double > xv( local_elems * nv );
312  success = NCFUNCAG( _vara_double )( _fileId, var_xv.varId, &startsv[0], &countsv[0], &xv[0] );
313  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for xv variable " );
314 
315  std::string yvstr( "yv" );
316  ReadNC::VarData& var_yv = _readNC->varInfo[yvstr];
317  std::vector< double > yv( local_elems * nv );
318  success = NCFUNCAG( _vara_double )( _fileId, var_yv.varId, &startsv[0], &countsv[0], &yv[0] );
319  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for yv variable " );
320 
321  // read other variables, like xc, yc, frac, area
322  std::string xcstr( "xc" );
323  ReadNC::VarData& var_xc = _readNC->varInfo[xcstr];
324  std::vector< double > xc( local_elems );
325  success = NCFUNCAG( _vara_double )( _fileId, var_xc.varId, &vmask.readStarts[0], &vmask.readCounts[0], &xc[0] );
326  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for xc variable " );
327 
328  std::string ycstr( "yc" );
329  ReadNC::VarData& var_yc = _readNC->varInfo[ycstr];
330  std::vector< double > yc( local_elems );
331  success = NCFUNCAG( _vara_double )( _fileId, var_yc.varId, &vmask.readStarts[0], &vmask.readCounts[0], &yc[0] );
332  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for yc variable " );
333 
334  std::string fracstr( "frac" );
335  ReadNC::VarData& var_frac = _readNC->varInfo[fracstr];
336  std::vector< double > frac( local_elems );
337  success = NCFUNCAG( _vara_double )( _fileId, var_frac.varId, &vmask.readStarts[0], &vmask.readCounts[0], &frac[0] );
338  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for frac variable " );
339  std::string areastr( "area" );
340  ReadNC::VarData& var_area = _readNC->varInfo[areastr];
341  std::vector< double > area( local_elems );
342  success = NCFUNCAG( _vara_double )( _fileId, var_area.varId, &vmask.readStarts[0], &vmask.readCounts[0], &area[0] );
343  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for area variable " );
344  // create tags for them
345  Tag areaTag, fracTag, xcTag, ycTag;
346  rval = mbImpl->tag_get_handle( "area", 1, MB_TYPE_DOUBLE, areaTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating area tag" );
347  rval = mbImpl->tag_get_handle( "frac", 1, MB_TYPE_DOUBLE, fracTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating frac tag" );
348  rval = mbImpl->tag_get_handle( "xc", 1, MB_TYPE_DOUBLE, xcTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating xc tag" );
349  rval = mbImpl->tag_get_handle( "yc", 1, MB_TYPE_DOUBLE, ycTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating yc tag" );
350 
351  //
352  EntityHandle* conn_arr;
353  EntityHandle vtx_handle;
354  Range tmp_range;
355 
356  // set connectivity into that space
357 
358  EntityHandle start_cell;
359  EntityType mdb_type = MBVERTEX;
360  if( nv == 3 )
361  mdb_type = MBTRI;
362  else if( nv == 4 )
363  mdb_type = MBQUAD;
364  else if( nv > 4 ) // (nv > 4)
365  mdb_type = MBPOLYGON;
366  // for nv = 1 , type is vertex
367 
368  if( nv > 1 && nb_with_mask1 > 0 )
369  {
370  rval = _readNC->readMeshIface->get_element_connect( nb_with_mask1, nv, mdb_type, 0, start_cell, conn_arr );MB_CHK_SET_ERR( rval, "Failed to create local cells" );
371  tmp_range.insert( start_cell, start_cell + nb_with_mask1 - 1 );
372  }
373 
374  // Create vertices; first identify different ones, with a tolerance
375  std::map< Node3D, EntityHandle > vertex_map;
376 
377  if( nb_with_mask1 > 0 )
378  {
379  // Set vertex coordinates
380  // will read all xv, yv, but use only those with correct mask on
381 
382  int elem_index = 0; // total index in netcdf arrays
383  const double pideg = acos( -1.0 ) / 180.0;
384 
385  for( ; elem_index < local_elems; elem_index++ )
386  {
387  if( 0 == mask[elem_index] ) continue; // nothing to do, do not advance elem_index in actual moab arrays
388  // set area and fraction on those elements too
389  for( int k = 0; k < nv; k++ )
390  {
391  int index_v_arr = nv * elem_index + k;
392  double x, y;
393  if( nv > 1 )
394  {
395  x = xv[index_v_arr];
396  y = yv[index_v_arr];
397  double cosphi = cos( pideg * y );
398  double zmult = sin( pideg * y );
399  double xmult = cosphi * cos( x * pideg );
400  double ymult = cosphi * sin( x * pideg );
401  Node3D pt( xmult, ymult, zmult );
402  vertex_map[pt] = 0;
403  }
404  else
405  {
406  x = xc[elem_index];
407  y = yc[elem_index];
408  Node3D pt( x, y, 0 );
409  vertex_map[pt] = 0;
410  }
411  }
412  }
413  int nLocalVertices = (int)vertex_map.size();
414  std::vector< double* > arrays;
415  EntityHandle start_vertex;
416  rval = _readNC->readMeshIface->get_node_coords( 3, nLocalVertices, 0, start_vertex, arrays );MB_CHK_SET_ERR( rval, "Failed to create local vertices" );
417 
418  vtx_handle = start_vertex;
419  // Copy vertex coordinates into entity sequence coordinate arrays
420  // and copy handle into vertex_map.
421  double *x = arrays[0], *y = arrays[1], *z = arrays[2];
422  for( auto i = vertex_map.begin(); i != vertex_map.end(); ++i )
423  {
424  i->second = vtx_handle;
425  ++vtx_handle;
426  *x = i->first.coords[0];
427  ++x;
428  *y = i->first.coords[1];
429  ++y;
430  *z = i->first.coords[2];
431  ++z;
432  }
433 
434  // int nj = gDims[4]-gDims[1]; // is it about 1 in irregular cases
435 
436  // int local_row_size = lCDims[3] - lCDims[0];
437  int global_row_size = gDims[3] - gDims[0]; // this is along
438  elem_index = -1;
439  int index = 0; // consider the mask for advancing in moab arrays;
440 
441  // create now vertex arrays, size vertex_map.size()
442  for( int j = lCDims[1]; j < lCDims[4]; j++ )
443  for( int i = lCDims[0]; i < lCDims[3]; i++ )
444  {
445  elem_index++;
446  if( 0 == mask[elem_index] ) continue; // nothing to do, do not advance elem_index in actual moab arrays
447  // set area and fraction on those elements too
448  for( int k = 0; k < nv; k++ )
449  {
450  int index_v_arr = nv * elem_index + k;
451  if( nv > 1 )
452  {
453  double x = xv[index_v_arr];
454  double y = yv[index_v_arr];
455  double cosphi = cos( pideg * y );
456  double zmult = sin( pideg * y );
457  double xmult = cosphi * cos( x * pideg );
458  double ymult = cosphi * sin( x * pideg );
459  Node3D pt( xmult, ymult, zmult );
460  conn_arr[index * nv + k] = vertex_map[pt];
461  }
462  }
463  EntityHandle cell = start_vertex + index;
464  if( nv > 1 ) cell = start_cell + index;
465  // set other tags, like xc, yc, frac, area
466  rval = mbImpl->tag_set_data( xcTag, &cell, 1, &xc[elem_index] );MB_CHK_SET_ERR( rval, "Failed to set xc tag" );
467  rval = mbImpl->tag_set_data( ycTag, &cell, 1, &yc[elem_index] );MB_CHK_SET_ERR( rval, "Failed to set yc tag" );
468  rval = mbImpl->tag_set_data( areaTag, &cell, 1, &area[elem_index] );MB_CHK_SET_ERR( rval, "Failed to set area tag" );
469  rval = mbImpl->tag_set_data( fracTag, &cell, 1, &frac[elem_index] );MB_CHK_SET_ERR( rval, "Failed to set frac tag" );
470 
471  // set the global id too:
472  int globalId = j * global_row_size + i + 1;
473  rval = mbImpl->tag_set_data( mGlobalIdTag, &cell, 1, &globalId );MB_CHK_SET_ERR( rval, "Failed to set global id tag" );
474  index++;
475  }
476 
477  rval = mbImpl->add_entities( _fileSet, tmp_range );MB_CHK_SET_ERR( rval, "Failed to add new cells to current file set" );
478 
479  // modify local file set, to merge coincident vertices, and to correct repeated vertices in elements
480  std::vector< Tag > tagList;
481  tagList.push_back( mGlobalIdTag );
482  tagList.push_back( xcTag );
483  tagList.push_back( ycTag );
484  tagList.push_back( areaTag );
485  tagList.push_back( fracTag );
486  rval = IntxUtils::remove_padded_vertices( mbImpl, _fileSet, tagList );MB_CHK_SET_ERR( rval, "Failed to remove duplicate vertices" );
487 
488  rval = mbImpl->get_entities_by_dimension( _fileSet, 2, faces );MB_CHK_ERR( rval );
489  Range all_verts;
490  rval = mbImpl->get_connectivity( faces, all_verts );MB_CHK_ERR( rval );
491  //printf(" range vert size :%ld \n", all_verts.size());
492  rval = mbImpl->add_entities( _fileSet, all_verts );MB_CHK_ERR( rval );
493 
494  // need to add adjacencies; TODO: fix this for all nc readers
495  // copy this logic from migrate mesh in par comm graph
496  Core* mb = (Core*)mbImpl;
497  AEntityFactory* adj_fact = mb->a_entity_factory();
498  if( !adj_fact->vert_elem_adjacencies() )
499  adj_fact->create_vert_elem_adjacencies();
500  else
501  {
502  for( Range::iterator it = faces.begin(); it != faces.end(); ++it )
503  {
504  EntityHandle eh = *it;
505  const EntityHandle* conn = NULL;
506  int num_nodes = 0;
507  rval = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
508  adj_fact->notify_create_entity( eh, conn, num_nodes );
509  }
510  }
511  }
512 
513 #ifdef MOAB_HAVE_MPI
514  ParallelComm*& myPcomm = _readNC->myPcomm;
515  if( myPcomm )
516  {
517  double tol = 1.e-12; // this is the same as static tolerance in NCHelper
518  ParallelMergeMesh pmm( myPcomm, tol );
519  rval = pmm.merge( _fileSet,
520  /* do not do local merge*/ false,
521  /* 2d cells*/ 2 );MB_CHK_SET_ERR( rval, "Failed to merge vertices in parallel" );
522 
523  // assign global ids only for vertices, cells have them fine
524  rval = myPcomm->assign_global_ids( _fileSet, /*dim*/ 0 );MB_CHK_ERR( rval );
525  }
526 #endif
527 
528  return MB_SUCCESS;
529 }
530 } // namespace moab