Mesh Oriented datABase  (version 5.6.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"
4 #include "moab/IntxMesh/IntxUtils.hpp"
5 #include "AEntityFactory.hpp"
6 #ifdef MOAB_HAVE_MPI
7 #include "moab/ParallelMergeMesh.hpp"
8 #endif
9 #ifdef MOAB_HAVE_ZOLTAN
10 #include "moab/ZoltanPartitioner.hpp"
11 #include "moab/TupleList.hpp"
12 #endif
13 
14 #include <cmath>
15 #include <sstream>
16 
17 namespace moab
18 {
19 
20 bool NCHelperDomain::can_read_file( ReadNC* readNC, int /*fileId*/ )
21 {
22  std::vector< std::string >& dimNames = readNC->dimNames;
23 
24  // If dimension names "n" AND "ni" AND "nj" AND "nv" exist then it should be the Domain grid
25  // some files do not have n, which should be just ni*nj
26  if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "ni" ) ) != dimNames.end() ) &&
27  ( std::find( dimNames.begin(), dimNames.end(), std::string( "nj" ) ) != dimNames.end() ) &&
28  ( std::find( dimNames.begin(), dimNames.end(), std::string( "nv" ) ) != dimNames.end() ) )
29  {
30  // we do not test anymore if it is an actual CAM file
31  return true;
32  }
33 
34  return false;
35 }
36 
37 ErrorCode NCHelperDomain::init_mesh_vals()
38 {
39  Interface*& mbImpl = _readNC->mbImpl;
40  std::vector< std::string >& dimNames = _readNC->dimNames;
41  std::vector< int >& dimLens = _readNC->dimLens;
42  std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
43  DebugOutput& dbgOut = _readNC->dbgOut;
44 #ifdef MOAB_HAVE_MPI
45  bool& isParallel = _readNC->isParallel;
46 #endif
47  int& partMethod = _readNC->partMethod;
48  ScdParData& parData = _readNC->parData;
49 
50  // Look for names of i/j dimensions
51  // First i
52  std::vector< std::string >::iterator vit;
53  unsigned int idx;
54  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ni" ) ) != dimNames.end() )
55  idx = vit - dimNames.begin();
56  else
57  {
58  MB_SET_ERR( MB_FAILURE, "Couldn't find 'ni' variable" );
59  }
60  iDim = idx;
61  gDims[0] = 0;
62  gDims[3] = dimLens[idx];
63 
64  // Then j
65  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "nj" ) ) != dimNames.end() )
66  idx = vit - dimNames.begin();
67  else
68  {
69  MB_SET_ERR( MB_FAILURE, "Couldn't find 'nj' variable" );
70  }
71  jDim = idx;
72  gDims[1] = 0;
73  gDims[4] = dimLens[idx]; // Add 2 for the pole points ? not needed
74 
75  // do not use gcdims ? or use only gcdims?
76 
77  // Try a truly 2D mesh
78  gDims[2] = -1;
79  gDims[5] = -1;
80 
81  // Get number of vertices per cell
82  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "nv" ) ) != dimNames.end() )
83  idx = vit - dimNames.begin();
84  else
85  {
86  MB_SET_ERR( MB_FAILURE, "Couldn't find 'nv' dimension" );
87  }
88  nvDim = idx;
89  nv = dimLens[idx];
90 
91  // Parse options to get subset
92  int rank = 0, procs = 1;
93 #ifdef MOAB_HAVE_MPI
94  if( isParallel )
95  {
96  ParallelComm*& myPcomm = _readNC->myPcomm;
97  rank = myPcomm->proc_config().proc_rank();
98  procs = myPcomm->proc_config().proc_size();
99  }
100 #endif
101  if( procs > 1 )
102  {
103  for( int i = 0; i < 6; i++ )
104  parData.gDims[i] = gDims[i];
105  parData.partMethod = partMethod;
106  int pdims[3];
107 
108  locallyPeriodic[0] = locallyPeriodic[1] = locallyPeriodic[2] = 0;
109  MB_CHK_ERR( ScdInterface::compute_partition( procs, rank, parData, lDims, locallyPeriodic, pdims ) );
110  for( int i = 0; i < 3; i++ )
111  parData.pDims[i] = pdims[i];
112 
113  dbgOut.tprintf( 1, "Partition: %dx%d (out of %dx%d)\n", lDims[3] - lDims[0], lDims[4] - lDims[1],
114  gDims[3] - gDims[0], gDims[4] - gDims[1] );
115  if( 0 == rank )
116  dbgOut.tprintf( 1, "Contiguous chunks of size %d bytes.\n",
117  8 * ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ) );
118  }
119  else
120  {
121  for( int i = 0; i < 6; i++ )
122  lDims[i] = gDims[i];
123  locallyPeriodic[0] = globallyPeriodic[0];
124  }
125 
126  // Now get actual coordinate values for vertices and cell centers
127  lCDims[0] = lDims[0];
128 
129  lCDims[3] = lDims[3];
130 
131  // For FV models, will always be non-periodic in j
132  lCDims[1] = lDims[1];
133  lCDims[4] = lDims[4];
134 
135 #if 0
136  // Resize vectors to store values later
137  if (-1 != lDims[0])
138  ilVals.resize(lDims[3] - lDims[0] + 1);
139  if (-1 != lCDims[0])
140  ilCVals.resize(lCDims[3] - lCDims[0] + 1);
141  if (-1 != lDims[1])
142  jlVals.resize(lDims[4] - lDims[1] + 1);
143  if (-1 != lCDims[1])
144  jlCVals.resize(lCDims[4] - lCDims[1] + 1);
145  if (nTimeSteps > 0)
146  tVals.resize(nTimeSteps);
147 #endif
148 
149  dbgOut.tprintf( 1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4] );
150  dbgOut.tprintf( 1, "%d elements, %d vertices\n", ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ),
151  ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ) * nv );
152 
153  // For each variable, determine the entity location type and number of levels
154  std::map< std::string, ReadNC::VarData >::iterator mit;
155  for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
156  {
157  ReadNC::VarData& vd = ( *mit ).second;
158 
159  // Default entLoc is ENTLOCSET
160  if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
161  {
162  if( ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) &&
163  ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) )
164  vd.entLoc = ReadNC::ENTLOCFACE;
165  else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jDim ) != vd.varDims.end() ) &&
166  ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) )
167  vd.entLoc = ReadNC::ENTLOCNSEDGE;
168  else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) &&
169  ( std::find( vd.varDims.begin(), vd.varDims.end(), iDim ) != vd.varDims.end() ) )
170  vd.entLoc = ReadNC::ENTLOCEWEDGE;
171  }
172 
173  // Default numLev is 0
174  if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels;
175  }
176 
177  std::vector< std::string > ijdimNames( 2 );
178  ijdimNames[0] = "__ni";
179  ijdimNames[1] = "__nj";
180  std::string tag_name;
181  Tag tagh;
182 
183  // __<dim_name>_LOC_MINMAX (for slon, slat, lon and lat)
184  for( unsigned int i = 0; i != ijdimNames.size(); i++ )
185  {
186  std::vector< int > val( 2, 0 );
187  if( ijdimNames[i] == "__ni" )
188  {
189  val[0] = lDims[0];
190  val[1] = lDims[3];
191  }
192  else if( ijdimNames[i] == "__nj" )
193  {
194  val[0] = lDims[1];
195  val[1] = lDims[4];
196  }
197 
198  std::stringstream ss_tag_name;
199  ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
200  tag_name = ss_tag_name.str();
201  MB_CHK_SET_ERR( mbImpl->tag_get_handle( tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh,
202  MB_TAG_SPARSE | MB_TAG_CREAT ),
203  "Trouble creating conventional tag " << tag_name );
204  MB_CHK_SET_ERR( mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] ),
205  "Trouble setting data to conventional tag " << tag_name );
206  dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
207  }
208 
209  // __<dim_name>_LOC_VALS (for slon, slat, lon and lat)
210  // Assume all have the same data type as lon (expected type is float or double)
211  switch( varInfo["xc"].varDataType )
212  {
213  case NC_FLOAT:
214  case NC_DOUBLE:
215  break;
216  default:
217  MB_SET_ERR( MB_FAILURE, "Unexpected variable data type for 'lon'" );
218  }
219 
220  // do not need conventional tags
221  Tag convTagsCreated = 0;
222  int def_val = 0;
223  MB_CHK_SET_ERR( mbImpl->tag_get_handle( "__CONV_TAGS_CREATED", 1, MB_TYPE_INTEGER, convTagsCreated,
224  MB_TAG_SPARSE | MB_TAG_CREAT, &def_val ),
225  "Trouble getting _CONV_TAGS_CREATED tag" );
226  int create_conv_tags_flag = 1;
227  MB_CHK_SET_ERR( mbImpl->tag_set_data( convTagsCreated, &_fileSet, 1, &create_conv_tags_flag ),
228  "Trouble setting _CONV_TAGS_CREATED tag" );
229 
230  return MB_SUCCESS;
231 }
232 
233 ErrorCode NCHelperDomain::create_mesh( Range& faces )
234 {
235  Interface*& mbImpl = _readNC->mbImpl;
236  // std::string& fileName = _readNC->fileName;
237  Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
238  // const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
239  DebugOutput& dbgOut = _readNC->dbgOut;
240 
241  bool& culling = _readNC->culling;
242  bool& repartition = _readNC->repartition;
243 
244  int success = 0;
245 
246  int local_elems = ( lDims[4] - lDims[1] ) * ( lDims[3] - lDims[0] );
247  dbgOut.tprintf( 1, "local cells: %d \n", local_elems );
248 
249  // count how many will be with mask 1 here
250  // basically, read the mask variable on the local elements;
251  std::string maskstr( "mask" );
252  ReadNC::VarData& vmask = _readNC->varInfo[maskstr];
253 
254  // mask is (nj, ni)
255  vmask.readStarts.push_back( lDims[1] );
256  vmask.readStarts.push_back( lDims[0] );
257  vmask.readCounts.push_back( lDims[4] - lDims[1] );
258  vmask.readCounts.push_back( lDims[3] - lDims[0] );
259  std::vector< int > mask( local_elems );
260  success = NCFUNCAG( _vara_int )( _fileId, vmask.varId, &vmask.readStarts[0], &vmask.readCounts[0], &mask[0] );
261  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read int data for mask variable " );
262 
263  std::vector< int > gids( local_elems );
264  int elem_index = 0;
265  int global_row_size = gDims[3] - gDims[0]; // this is along first dimension in global decomposition
266  // create global id array for cells, for all cells, including those with 0 mask; which will be not used eventually
267  for( int j = lCDims[1]; j < lCDims[4]; j++ )
268  for( int i = lCDims[0]; i < lCDims[3]; i++ )
269  {
270  gids[elem_index] = j * global_row_size + i + 1;
271  elem_index++;
272  }
273  std::vector< NCDF_SIZE > startsv( 3 );
274  startsv[0] = vmask.readStarts[0];
275  startsv[1] = vmask.readStarts[1];
276  startsv[2] = 0;
277  std::vector< NCDF_SIZE > countsv( 3 );
278  countsv[0] = vmask.readCounts[0];
279  countsv[1] = vmask.readCounts[1];
280  countsv[2] = nv; // number of vertices per element
281 
282  // read xv and yv coords for vertices, and create elements;
283  std::string xvstr( "xv" );
284  ReadNC::VarData& var_xv = _readNC->varInfo[xvstr];
285  std::vector< double > xv( local_elems * nv );
286 
287  std::vector< NCDF_SIZE > starts_vv, counts_vv;
288  starts_vv.resize( 3 );
289  counts_vv.resize( 3 );
290  bool nv_last = true;
291  if( _readNC->dimNames[var_xv.varDims[0]] == std::string( "nv" ) ) // it means that nv is the first dimension
292  {
293  starts_vv[0] = 0;
294  starts_vv[1] = startsv[0];
295  starts_vv[2] = startsv[1];
296  counts_vv[0] = nv;
297  counts_vv[1] = countsv[0];
298  counts_vv[2] = countsv[1];
299  nv_last = false;
300  }
301  else
302  {
303  starts_vv = startsv;
304  counts_vv = countsv;
305  }
306  dbgOut.tprintf( 1, " nv is the last dimension in xv array (0 or 1) %d \n", (int)nv_last );
307  success = NCFUNCAG( _vara_double )( _fileId, var_xv.varId, &starts_vv[0], &counts_vv[0], &xv[0] );
308  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for xv variable " );
309 
310  std::string yvstr( "yv" );
311  ReadNC::VarData& var_yv = _readNC->varInfo[yvstr];
312  std::vector< double > yv( local_elems * nv );
313  success = NCFUNCAG( _vara_double )( _fileId, var_yv.varId, &starts_vv[0], &counts_vv[0], &yv[0] );
314  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for yv variable " );
315 
316  // read other variables, like xc, yc, frac, area
317  std::string xcstr( "xc" );
318  ReadNC::VarData& var_xc = _readNC->varInfo[xcstr];
319  std::vector< double > xc( local_elems );
320  success = NCFUNCAG( _vara_double )( _fileId, var_xc.varId, &vmask.readStarts[0], &vmask.readCounts[0], &xc[0] );
321  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for xc variable " );
322 
323  std::string ycstr( "yc" );
324  ReadNC::VarData& var_yc = _readNC->varInfo[ycstr];
325  std::vector< double > yc( local_elems );
326  success = NCFUNCAG( _vara_double )( _fileId, var_yc.varId, &vmask.readStarts[0], &vmask.readCounts[0], &yc[0] );
327  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for yc variable " );
328 
329  std::string fracstr( "frac" );
330  ReadNC::VarData& var_frac = _readNC->varInfo[fracstr];
331  std::vector< double > frac( local_elems );
332  if( var_frac.varId >= 0 )
333  {
334  success =
335  NCFUNCAG( _vara_double )( _fileId, var_frac.varId, &vmask.readStarts[0], &vmask.readCounts[0], &frac[0] );
336  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for frac variable " );
337  }
338  std::string areastr( "area" );
339  std::vector< double > area( local_elems );
340  ReadNC::VarData& var_area = _readNC->varInfo[areastr];
341  success = NCFUNCAG( _vara_double )( _fileId, var_area.varId, &vmask.readStarts[0], &vmask.readCounts[0], &area[0] );
342  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for area variable " );
343  // create tags for them
344  Tag areaTag, fracTag, xcTag, ycTag;
345  MB_CHK_SET_ERR( mbImpl->tag_get_handle( "area", 1, MB_TYPE_DOUBLE, areaTag, MB_TAG_DENSE | MB_TAG_CREAT ),
346  "Trouble creating area tag" );
347  MB_CHK_SET_ERR( mbImpl->tag_get_handle( "frac", 1, MB_TYPE_DOUBLE, fracTag, MB_TAG_DENSE | MB_TAG_CREAT ),
348  "Trouble creating frac tag" );
349  MB_CHK_SET_ERR( mbImpl->tag_get_handle( "xc", 1, MB_TYPE_DOUBLE, xcTag, MB_TAG_DENSE | MB_TAG_CREAT ),
350  "Trouble creating xc tag" );
351  MB_CHK_SET_ERR( mbImpl->tag_get_handle( "yc", 1, MB_TYPE_DOUBLE, ycTag, MB_TAG_DENSE | MB_TAG_CREAT ),
352  "Trouble creating yc tag" );
353 
354  // create tags for GRID_IMASK, which will be the same name as the Scrip helper tag that holds the mask
355  // create the maskTag GRID_IMASK, with default value of 1
356  Tag maskTag;
357  int def_val = 1;
358  MB_CHK_SET_ERR( mbImpl->tag_get_handle( "GRID_IMASK", 1, MB_TYPE_INTEGER, maskTag, MB_TAG_DENSE | MB_TAG_CREAT,
359  &def_val ),
360  "Trouble creating GRID_IMASK tag" );
361 
362  // will now look to repartition the cells, using zoltan, looking at the xc and yc coordinates in 2d, convert to 3d,
363  // and decide based on those partitioning info to what task to send each cell, along with its vertices, and global id
364 #ifdef MOAB_HAVE_MPI
365  int procs = 1;
366  bool& isParallel = _readNC->isParallel;
367  ParallelComm* myPcomm = NULL;
368  if( isParallel )
369  {
370  myPcomm = _readNC->myPcomm;
371  procs = myPcomm->proc_config().proc_size();
372  }
373 
374  if( procs >= 2 && repartition )
375  {
376  // Redistribute local cells after trivial partition (e.g. apply Zoltan partition)
377  MB_CHK_SET_ERR( redistribute_cells( myPcomm, xc, yc, xv, yv, frac, mask, area, gids, nv, nv_last ),
378  "Failed to redistribute local cells" );
379  local_elems = (int)xc.size();
380  dbgOut.tprintf( 1, "local cells after repartition: %d \n", local_elems );
381  }
382 
383 #endif
384 
385  int nb_with_mask1 = 0;
386  for( int i = 0; i < local_elems; i++ )
387  if( 1 == mask[i] ) nb_with_mask1++;
388  dbgOut.tprintf( 1, "local cells with mask 1: %d \n", nb_with_mask1 );
389 
390  EntityHandle* conn_arr;
391  EntityHandle vtx_handle;
392  Range tmp_range;
393 
394  // set connectivity into that space
395 
396  EntityHandle start_cell;
397  EntityType mdb_type = MBVERTEX;
398  if( nv == 3 )
399  mdb_type = MBTRI;
400  else if( nv == 4 )
401  mdb_type = MBQUAD;
402  else if( nv > 4 ) // (nv > 4)
403  mdb_type = MBPOLYGON;
404  // for nv = 1 , type is vertex
405 
406  int num_actual_cells = local_elems;
407  if( culling ) num_actual_cells = nb_with_mask1;
408 
409  if( nv > 1 && num_actual_cells > 0 )
410  {
411  MB_CHK_SET_ERR( _readNC->readMeshIface->get_element_connect( num_actual_cells, nv, mdb_type, 0, start_cell,
412  conn_arr ),
413  "Failed to create local cells" );
414  tmp_range.insert( start_cell, start_cell + num_actual_cells - 1 );
415  }
416 
417  // Create vertices; first identify different ones, with a tolerance
418  std::map< Node3D, EntityHandle > vertex_map;
419 
420  if( num_actual_cells > 0 )
421  {
422  // Set vertex coordinates
423  // will read all xv, yv, but use only those with correct mask on
424 
425  // total index in netcdf arrays
426  const double pideg = acos( -1.0 ) / 180.0;
427 
428  for( elem_index = 0; elem_index < local_elems; elem_index++ )
429  {
430  if( culling && 0 == mask[elem_index] )
431  continue; // nothing to do, do not advance elem_index in actual moab arrays
432  // set area and fraction on those elements too
433  dbgOut.tprintf( 3, "elem index %d \n", elem_index );
434  for( int k = 0; k < nv; k++ )
435  {
436  int index_v_arr = nv * elem_index + k;
437  if( !nv_last ) index_v_arr = k * local_elems + elem_index;
438  double x, y;
439  if( nv > 1 )
440  {
441  x = xv[index_v_arr];
442  y = yv[index_v_arr];
443  double cosphi = cos( pideg * y );
444  double zmult = sin( pideg * y );
445  double xmult = cosphi * cos( x * pideg );
446  double ymult = cosphi * sin( x * pideg );
447  Node3D pt( xmult, ymult, zmult );
448  vertex_map[pt] = 0;
449  dbgOut.tprintf( 3, " %d x=%5.1f y=%5.1f pt: %f %f %f \n", k, x, y, pt.coords[0], pt.coords[1],
450  pt.coords[2] );
451  }
452  else
453  {
454  x = xc[elem_index];
455  y = yc[elem_index];
456  Node3D pt( x, y, 0 );
457  vertex_map[pt] = 0;
458  }
459  }
460  }
461  int nLocalVertices = (int)vertex_map.size();
462  std::vector< double* > arrays;
463  EntityHandle start_vertex;
464  MB_CHK_SET_ERR( _readNC->readMeshIface->get_node_coords( 3, nLocalVertices, 0, start_vertex, arrays ),
465  "Failed to create local vertices" );
466 
467  vtx_handle = start_vertex;
468  // Copy vertex coordinates into entity sequence coordinate arrays
469  // and copy handle into vertex_map.
470  double *x = arrays[0], *y = arrays[1], *z = arrays[2];
471  for( auto i = vertex_map.begin(); i != vertex_map.end(); ++i )
472  {
473  i->second = vtx_handle;
474  ++vtx_handle;
475  *x = i->first.coords[0];
476  ++x;
477  *y = i->first.coords[1];
478  ++y;
479  *z = i->first.coords[2];
480  ++z;
481  }
482 
483  // int nj = gDims[4]-gDims[1]; // is it about 1 in irregular cases
484 
485  // int local_row_size = lCDims[3] - lCDims[0];
486  int index = 0; // consider the mask for advancing in moab arrays;
487 
488  // create now vertex arrays, size vertex_map.size()
489  for( elem_index = 0; elem_index < local_elems; elem_index++ )
490  {
491  if( culling && 0 == mask[elem_index] )
492  continue; // nothing to do, do not advance elem_index in actual moab arrays
493  // set area and fraction on those elements too
494  for( int k = 0; k < nv; k++ )
495  {
496  int index_v_arr = nv * elem_index + k;
497  if( !nv_last ) index_v_arr = k * local_elems + elem_index;
498  if( nv > 1 )
499  {
500  double x = xv[index_v_arr];
501  double y = yv[index_v_arr];
502  double cosphi = cos( pideg * y );
503  double zmult = sin( pideg * y );
504  double xmult = cosphi * cos( x * pideg );
505  double ymult = cosphi * sin( x * pideg );
506  Node3D pt( xmult, ymult, zmult );
507  conn_arr[index * nv + k] = vertex_map[pt];
508  }
509  }
510  EntityHandle cell = start_vertex + index;
511  if( nv > 1 ) cell = start_cell + index;
512  // set other tags, like xc, yc, frac, area
513  MB_CHK_SET_ERR( mbImpl->tag_set_data( xcTag, &cell, 1, &xc[elem_index] ), "Failed to set xc tag" );
514  MB_CHK_SET_ERR( mbImpl->tag_set_data( ycTag, &cell, 1, &yc[elem_index] ), "Failed to set yc tag" );
515  MB_CHK_SET_ERR( mbImpl->tag_set_data( areaTag, &cell, 1, &area[elem_index] ), "Failed to set area tag" );
516  MB_CHK_SET_ERR( mbImpl->tag_set_data( fracTag, &cell, 1, &frac[elem_index] ), "Failed to set frac tag" );
517  MB_CHK_SET_ERR( mbImpl->tag_set_data( maskTag, &cell, 1, &mask[elem_index] ), "Failed to set mask tag" );
518 
519  // set the global id too:
520  int globalId = gids[elem_index];
521  MB_CHK_SET_ERR( mbImpl->tag_set_data( mGlobalIdTag, &cell, 1, &globalId ), "Failed to set global id tag" );
522  index++;
523  }
524 
525  MB_CHK_SET_ERR( mbImpl->add_entities( _fileSet, tmp_range ), "Failed to add new cells to current file set" );
526 
527  // modify local file set, to merge coincident vertices, and to correct repeated vertices in elements
528  std::vector< Tag > tagList;
529  tagList.push_back( mGlobalIdTag );
530  tagList.push_back( xcTag );
531  tagList.push_back( ycTag );
532  tagList.push_back( areaTag );
533  tagList.push_back( fracTag );
534  tagList.push_back( maskTag ); // not sure this is needed though ? on cells or on vertices?
535  MB_CHK_SET_ERR( IntxUtils::remove_padded_vertices( mbImpl, _fileSet, tagList ),
536  "Failed to remove duplicate vertices" );
537 
538  MB_CHK_ERR( mbImpl->get_entities_by_dimension( _fileSet, 2, faces ) );
539  Range all_verts;
540  MB_CHK_ERR( mbImpl->get_connectivity( faces, all_verts ) );
541  //printf(" range vert size :%ld \n", all_verts.size());
542  MB_CHK_ERR( mbImpl->add_entities( _fileSet, all_verts ) );
543 
544  // need to add adjacencies; TODO: fix this for all nc readers
545  // copy this logic from migrate mesh in par comm graph
546  Core* mb = (Core*)mbImpl;
547  AEntityFactory* adj_fact = mb->a_entity_factory();
548  if( !adj_fact->vert_elem_adjacencies() )
549  adj_fact->create_vert_elem_adjacencies();
550  else
551  {
552  for( Range::iterator it = faces.begin(); it != faces.end(); ++it )
553  {
554  EntityHandle eh = *it;
555  const EntityHandle* conn = NULL;
556  int num_nodes = 0;
557  MB_CHK_ERR( mb->get_connectivity( eh, conn, num_nodes ) );
558  adj_fact->notify_create_entity( eh, conn, num_nodes );
559  }
560  }
561  }
562 
563 #ifdef MOAB_HAVE_MPI
564  myPcomm = _readNC->myPcomm; // we will have to set the global id on vertices in any case, even
565  if( myPcomm && procs >= 2 )
566  {
567  double tol = 1.e-12; // this is the same as static tolerance in NCHelper
568  ParallelMergeMesh pmm( myPcomm, tol );
569  MB_CHK_SET_ERR( pmm.merge( _fileSet,
570  /* do not do local merge*/ false,
571  /* 2d cells*/ 2 ),
572  "Failed to merge vertices in parallel" );
573  // assign global ids only for vertices, cells have them fine
574  MB_CHK_ERR( myPcomm->assign_global_ids( _fileSet, /*dim*/ 0 ) );
575  }
576 #endif
577 
578  return MB_SUCCESS;
579 }
580 
581 #ifdef MOAB_HAVE_MPI
582 ErrorCode NCHelperDomain::redistribute_cells( ParallelComm* myPcomm,
583  std::vector< double >& xc, // center x
584  std::vector< double >& yc, // center y
585  std::vector< double >& xv, // vertex coords
586  std::vector< double >& yv,
587  std::vector< double >& frac, // fractions
588  std::vector< int >& masks, // mask
589  std::vector< double >& area, // area
590  std::vector< int >& gids, // global ids
591  int nv, // number of vertices per cell
592  bool nv_last ) // type of xv, yv, first or last
593 {
594  const int my_rank = myPcomm->proc_config().proc_rank();
595 
596 #ifdef MOAB_HAVE_ZOLTAN
597 
598  if( !my_rank ) std::cout << "Using Zoltan RCB spatial partitioning method\n";
599 
600  size_t num_local_cells = gids.size();
601  if( _readNC->culling )
602  {
603  num_local_cells = std::count( masks.begin(), masks.end(), 1 );
604  }
605 
606  size_t actual_index = 0;
607  std::vector< double > xi( num_local_cells ), yi( num_local_cells ), zi( num_local_cells );
608  std::vector< int > gids2( num_local_cells );
609 
610  // now loop over coordinates and accumulate
611  const double deg_to_rad = std::acos( -1.0 ) / 180.0;
612  for( size_t i = 0; i < xc.size(); ++i )
613  {
614  if( _readNC->culling && masks[i] == 0 ) continue;
615 
616  const double lon_rad = xc[i] * deg_to_rad;
617  const double lat_rad = yc[i] * deg_to_rad;
618  const double cos_lat = std::cos( lat_rad );
619  const double sin_lat = std::sin( lat_rad );
620 
621  xi[actual_index] = cos_lat * std::cos( lon_rad );
622  yi[actual_index] = cos_lat * std::sin( lon_rad );
623  zi[actual_index] = sin_lat;
624  gids2[actual_index] = gids[i];
625 
626  ++actual_index;
627  }
628 
629  std::vector< int > dest( num_local_cells );
630  {
631  // apply Zoltan partitioning using RCB strategy
632  ZoltanPartitioner mbZTool( _readNC->mbImpl, myPcomm, false, 0, NULL );
633  MB_CHK_SET_ERR( mbZTool.repartition_to_procs( xi, yi, zi, gids2, "RCB", dest ),
634  "Error in Zoltan partitioning" );
635  }
636 
637  // now use crystal router to send the arrays to the right places
638  moab::TupleList tl;
639  const unsigned num_real = 2 * nv + 4;
640  tl.initialize( 3, 0, 0, num_real, num_local_cells ); // to_proc, global_id, mask
641  tl.enableWriteAccess();
642 
643  // populate
644  size_t index_in_dest = 0;
645  for( size_t i = 0; i < xc.size(); ++i )
646  {
647  if( _readNC->culling && masks[i] == 0 ) continue;
648 
649  const int to_proc = dest[index_in_dest];
650  const int global_id = gids2[index_in_dest];
651  const int mask = masks[i];
652  const int n = tl.get_n();
653 
654  tl.vi_wr[3 * n + 0] = to_proc;
655  tl.vi_wr[3 * n + 1] = global_id;
656  tl.vi_wr[3 * n + 2] = mask;
657 
658  tl.vr_wr[n * num_real + 0] = area[i];
659  tl.vr_wr[n * num_real + 1] = xc[i];
660  tl.vr_wr[n * num_real + 2] = yc[i];
661  tl.vr_wr[n * num_real + 3] = frac[i];
662 
663  for( int k = 0; k < nv; ++k )
664  {
665  const size_t index_v = nv_last ? nv * i + k : k * xc.size() + i;
666  tl.vr_wr[n * num_real + 4 + k] = xv[index_v];
667  tl.vr_wr[n * num_real + 4 + nv + k] = yv[index_v];
668  }
669 
670  tl.inc_n();
671  ++index_in_dest;
672  }
673 
674  // Communicate using crystal router
675  myPcomm->proc_config().crystal_router()->gs_transfer( 1, tl, 0 );
676 
677  // after communication, on each processor we should have tuples coming in
678  // rearrange (sort) the vectors by global id
679  moab::TupleList::buffer sort_buffer;
680  int N = tl.get_n();
681  sort_buffer.buffer_init( N );
682  tl.sort( 1, &sort_buffer ); // sort by global ID (index 1)
683 
684  // Resize and fill output arrays
685  xc.resize( N );
686  yc.resize( N );
687  xv.resize( N * nv );
688  yv.resize( N * nv );
689  frac.resize( N );
690  masks.resize( N );
691  area.resize( N );
692  gids.resize( N );
693 
694  for( int n = 0; n < N; ++n )
695  {
696  gids[n] = tl.vi_wr[3 * n + 1];
697  masks[n] = tl.vi_wr[3 * n + 2];
698  area[n] = tl.vr_wr[n * num_real + 0];
699  xc[n] = tl.vr_wr[n * num_real + 1];
700  yc[n] = tl.vr_wr[n * num_real + 2];
701  frac[n] = tl.vr_wr[n * num_real + 3];
702 
703  for( int k = 0; k < nv; ++k )
704  {
705  const size_t index_v = nv_last ? nv * n + k : k * N + n;
706  xv[index_v] = tl.vr_wr[n * num_real + 4 + k];
707  yv[index_v] = tl.vr_wr[n * num_real + 4 + nv + k];
708  }
709  }
710 
711  return MB_SUCCESS;
712 #else
713  const int num_procs = myPcomm->proc_config().proc_size();
714  // Trivial partitioning fallback if Zoltan is not available
715  constexpr bool use_spatial_partitioning = true; // <-- Toggle this flag to test spatial partitioning
716 
717  if( !my_rank )
718  {
719  if( use_spatial_partitioning && !_readNC->culling )
720  std::cout << "Using native spatial partitioning method\n";
721  else
722  std::cout << "Using native trivial partitioning method\n";
723  }
724 
725  size_t num_local_cells = gids.size();
726  if( _readNC->culling )
727  {
728  num_local_cells = std::count( masks.begin(), masks.end(), 1 );
729  }
730 
731  moab::TupleList tl;
732  const unsigned num_real = 2 * nv + 4;
733  tl.initialize( 3, 0, 0, num_real, num_local_cells ); // to_proc, global_id, mask
734  tl.enableWriteAccess();
735 
736  auto clamp = []( int v, int lo, int hi ) { return ( v < lo ) ? lo : ( v > hi ) ? hi : v; };
737 
738  size_t actual_index = 0;
739  for( size_t i = 0; i < xc.size(); ++i )
740  {
741  if( _readNC->culling && masks[i] == 0 ) continue;
742 
743  const int n = tl.get_n();
744  int to_proc = 0;
745  if( use_spatial_partitioning && !_readNC->culling ) // probably will not work when culling is on
746  {
747  // Spatial partitioning using longitude in range [-180, 180]
748  double lon = xc[i];
749  // Normalize longitude to [-180, 180]
750  while( lon < -180.0 )
751  lon += 360.0;
752  while( lon > 180.0 )
753  lon -= 360.0;
754  // Now map to [0, 1]
755  double lon_norm = ( lon + 180.0 ) / 360.0;
756 
757  // Now let us compute the right partition
758  to_proc = static_cast< int >( lon_norm * num_procs );
759  // Handle edge case where lon_norm == 1.0 (e.g. xc[i] == 180)
760  to_proc = clamp( to_proc, 0, num_procs - 1 );
761  }
762  else
763  {
764  // Trivial linear partitioner
765  // to_proc = gids[i] % num_procs;
766  // Block-based partitioning: roughly equal number of cells per processor
767  to_proc = actual_index * num_procs / num_local_cells;
768  }
769 
770  tl.vi_wr[3 * n + 0] = to_proc;
771  tl.vi_wr[3 * n + 1] = gids[i];
772  tl.vi_wr[3 * n + 2] = masks[i];
773 
774  tl.vr_wr[n * num_real + 0] = area[i];
775  tl.vr_wr[n * num_real + 1] = xc[i];
776  tl.vr_wr[n * num_real + 2] = yc[i];
777  tl.vr_wr[n * num_real + 3] = frac[i];
778 
779  for( int k = 0; k < nv; ++k )
780  {
781  const size_t index_v = nv_last ? nv * i + k : k * xc.size() + i;
782  tl.vr_wr[n * num_real + 4 + k] = xv[index_v];
783  tl.vr_wr[n * num_real + 4 + nv + k] = yv[index_v];
784  }
785 
786  tl.inc_n();
787  ++actual_index;
788  }
789 
790  // Communicate using crystal router
791  myPcomm->proc_config().crystal_router()->gs_transfer( 1, tl, 0 );
792 
793  // Sort by global ID
794  moab::TupleList::buffer sort_buffer;
795  const int N = tl.get_n();
796  sort_buffer.buffer_init( N );
797  tl.sort( 1, &sort_buffer ); // sort by global ID
798 
799  // Resize and fill output arrays
800  xc.resize( N );
801  yc.resize( N );
802  xv.resize( N * nv );
803  yv.resize( N * nv );
804  frac.resize( N );
805  masks.resize( N );
806  area.resize( N );
807  gids.resize( N );
808 
809  for( int n = 0; n < N; ++n )
810  {
811  gids[n] = tl.vi_wr[3 * n + 1];
812  masks[n] = tl.vi_wr[3 * n + 2];
813  area[n] = tl.vr_wr[n * num_real + 0];
814  xc[n] = tl.vr_wr[n * num_real + 1];
815  yc[n] = tl.vr_wr[n * num_real + 2];
816  frac[n] = tl.vr_wr[n * num_real + 3];
817 
818  for( int k = 0; k < nv; ++k )
819  {
820  const size_t index_v = nv_last ? nv * n + k : k * N + n;
821  xv[index_v] = tl.vr_wr[n * num_real + 4 + k];
822  yv[index_v] = tl.vr_wr[n * num_real + 4 + nv + k];
823  }
824  }
825 
826  return MB_SUCCESS;
827 #endif
828 }
829 
830 #endif
831 
832 } // namespace moab