 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /*
2 * NCHelperESMF.cpp
3 *
4 * Created on: Sep. 12, 2023
5 */
6
7 #include "NCHelperESMF.hpp"
8 #include "moab/ReadUtilIface.hpp"
9 #include "moab/FileOptions.hpp"
10 #include "MBTagConventions.hpp"
11
12 #ifdef MOAB_HAVE_ZOLTAN
13 #include "moab/ZoltanPartitioner.hpp"
14 #endif
15
16 #include <cmath>
17
18 namespace moab
19 {
20
21 const int DEFAULT_MAX_EDGES_PER_CELL = 10;
22 const double pideg = acos( -1.0 ) / 180.0;
23
24 NCHelperESMF::NCHelperESMF( ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet )
25 : UcdNCHelper( readNC, fileId, opts, fileSet ), maxEdgesPerCell( DEFAULT_MAX_EDGES_PER_CELL ), coordDim( 0 ),
26 centerCoordsId( -1 ), degrees( true )
27 {
28 }
29
30 bool NCHelperESMF::can_read_file( ReadNC* readNC )
31 {
32 std::vector< std::string >& dimNames = readNC->dimNames;
33 if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "nodeCount" ) ) != dimNames.end() ) &&
34 ( std::find( dimNames.begin(), dimNames.end(), std::string( "elementCount" ) ) != dimNames.end() ) &&
35 ( std::find( dimNames.begin(), dimNames.end(), std::string( "maxNodePElement" ) ) != dimNames.end() ) &&
36 ( std::find( dimNames.begin(), dimNames.end(), std::string( "coordDim" ) ) != dimNames.end() ) )
37 {
38 return true;
39 }
40
41 return false;
42 }
43
44 ErrorCode NCHelperESMF::init_mesh_vals()
45 {
46 std::vector< std::string >& dimNames = _readNC->dimNames;
47 std::vector< int >& dimLens = _readNC->dimLens;
48 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
49
50 ErrorCode rval;
51 unsigned int idx;
52 std::vector< std::string >::iterator vit;
53
54 // Get max edges per cell reported in the ESMF file header
55 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "maxNodePElement" ) ) != dimNames.end() )
56 {
57 idx = vit - dimNames.begin();
58 maxEdgesPerCell = dimLens[idx];
59 if( maxEdgesPerCell > DEFAULT_MAX_EDGES_PER_CELL )
60 {
61 MB_SET_ERR( MB_INVALID_SIZE,
62 "maxEdgesPerCell read from the ESMF file header has exceeded " << DEFAULT_MAX_EDGES_PER_CELL );
63 }
64 }
65
66 // Get number of cells
67 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "elementCount" ) ) != dimNames.end() )
68 idx = vit - dimNames.begin();
69 else
70 {
71 MB_SET_ERR( MB_FAILURE, "Couldn't find 'elementCount' dimension" );
72 }
73 cDim = idx;
74 nCells = dimLens[idx];
75
76 // Get number of vertices per cell (max)
77 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "maxNodePElement" ) ) != dimNames.end() )
78 idx = vit - dimNames.begin();
79 else
80 {
81 MB_SET_ERR( MB_FAILURE, "Couldn't find 'maxNodePElement' dimension" );
82 }
83
84 maxEdgesPerCell = dimLens[idx];
85
86 // Get number of vertices
87 vDim = -1;
88 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "nodeCount" ) ) != dimNames.end() )
89 idx = vit - dimNames.begin();
90 else
91 {
92 MB_SET_ERR( MB_FAILURE, "Couldn't find 'nodeCount' dimension" );
93 }
94 vDim = idx;
95 nVertices = dimLens[idx];
96
97 // Get number of coord dimensions
98 if( ( vit = std::find( dimNames.begin(), dimNames.end(), "coordDim" ) ) != dimNames.end() )
99 idx = vit - dimNames.begin();
100 else
101 {
102 MB_SET_ERR( MB_FAILURE, "Couldn't find 'coordDim' dimension" );
103 }
104
105 coordDim = dimLens[idx];
106
107 int success = NCFUNC( inq_varid )( _fileId, "centerCoords", ¢erCoordsId );
108 if( success ) centerCoordsId = -1; // no center coords variable
109
110 // decide now the units, by looking at nodeCoords; they should always exist
111 int nodeCoordsId;
112 success = NCFUNC( inq_varid )( _fileId, "nodeCoords", &nodeCoordsId );
113 if( success ) MB_CHK_SET_ERR( MB_FAILURE, "Trouble getting nodeCoords" );
114
115 auto vmit = varInfo.find( "nodeCoords" );
116 if( varInfo.end() == vmit )
117 MB_SET_ERR( MB_FAILURE, "Couldn't find variable "
118 << "nodeCoords" );
119 ReadNC::VarData& glData = vmit->second;
120 auto attIt = glData.varAtts.find( "units" );
121 if( attIt != glData.varAtts.end() )
122 {
123 unsigned int sz = attIt->second.attLen;
124 std::string att_data;
125 att_data.resize( sz + 1 );
126 att_data[sz] = '\000';
127 success =
128 NCFUNC( get_att_text )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
129 if( 0 == success && att_data.find( "radians" ) != std::string::npos ) degrees = false;
130 }
131
132 // Hack: create dummy variables for dimensions (like nCells) with no corresponding coordinate
133 // variables
134 rval = create_dummy_variables();MB_CHK_SET_ERR( rval, "Failed to create dummy variables" );
135
136 return MB_SUCCESS;
137 }
138
139 ErrorCode NCHelperESMF::create_mesh( Range& faces )
140 {
141 bool& noMixedElements = _readNC->noMixedElements;
142 DebugOutput& dbgOut = _readNC->dbgOut;
143
144 #ifdef MOAB_HAVE_MPI
145 int rank = 0;
146 int procs = 1;
147 bool& isParallel = _readNC->isParallel;
148 ParallelComm* myPcomm = NULL;
149 if( isParallel )
150 {
151 myPcomm = _readNC->myPcomm;
152 rank = myPcomm->proc_config().proc_rank();
153 procs = myPcomm->proc_config().proc_size();
154 }
155
156 if( procs >= 2 )
157 {
158 // Shift rank to obtain a rotated trivial partition
159 int shifted_rank = rank;
160 int& trivialPartitionShift = _readNC->trivialPartitionShift;
161 if( trivialPartitionShift > 0 ) shifted_rank = ( rank + trivialPartitionShift ) % procs;
162
163 // Compute the number of local cells on this proc
164 nLocalCells = int( std::floor( 1.0 * nCells / procs ) );
165
166 // The starting global cell index in the ESMF file for this proc
167 int start_cell_idx = shifted_rank * nLocalCells;
168
169 // Number of extra cells after equal split over procs
170 int iextra = nCells % procs;
171
172 // Allocate extra cells over procs
173 if( shifted_rank < iextra ) nLocalCells++;
174 start_cell_idx += std::min( shifted_rank, iextra );
175
176 start_cell_idx++; // 0 based -> 1 based
177
178 // Redistribute local cells after trivial partition (e.g. apply Zoltan partition)
179 ErrorCode rval = redistribute_local_cells( start_cell_idx, myPcomm );MB_CHK_SET_ERR( rval, "Failed to redistribute local cells after trivial partition" );
180 }
181 else
182 {
183 nLocalCells = nCells;
184 localGidCells.insert( 1, nLocalCells );
185 }
186 #else
187 nLocalCells = nCells;
188 localGidCells.insert( 1, nLocalCells );
189 #endif
190 dbgOut.tprintf( 1, " localGidCells.psize() = %d\n", (int)localGidCells.psize() );
191 dbgOut.tprintf( 1, " localGidCells.size() = %d\n", (int)localGidCells.size() );
192
193 // Read number of edges on each local cell, to calculate actual maxEdgesPerCell
194 int nEdgesOnCellVarId;
195 int success = NCFUNC( inq_varid )( _fileId, "numElementConn", &nEdgesOnCellVarId );
196 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of numElementConn" );
197 std::vector< int > num_edges_on_local_cells( nLocalCells );
198
199 #ifdef MOAB_HAVE_PNETCDF
200 size_t nb_reads = localGidCells.psize();
201 std::vector< int > requests( nb_reads );
202 std::vector< int > statuss( nb_reads );
203 size_t idxReq = 0;
204 #endif
205 size_t indexInArray = 0;
206 for( Range::pair_iterator pair_iter = localGidCells.pair_begin(); pair_iter != localGidCells.pair_end();
207 ++pair_iter )
208 {
209 EntityHandle starth = pair_iter->first;
210 EntityHandle endh = pair_iter->second;
211 NCDF_SIZE read_start = (NCDF_SIZE)( starth - 1 );
212 NCDF_SIZE read_count = (NCDF_SIZE)( endh - starth + 1 );
213
214 // Do a partial read in each subrange
215 #ifdef MOAB_HAVE_PNETCDF
216 success = NCFUNCREQG( _vara_int )( _fileId, nEdgesOnCellVarId, &read_start, &read_count,
217 &( num_edges_on_local_cells[indexInArray] ), &requests[idxReq++] );
218 #else
219 success = NCFUNCAG( _vara_int )( _fileId, nEdgesOnCellVarId, &read_start, &read_count,
220 &( num_edges_on_local_cells[indexInArray] ) );
221 #endif
222 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read nEdgesOnCell data in a loop" );
223
224 // Increment the index for next subrange
225 indexInArray += ( endh - starth + 1 );
226 }
227
228 #ifdef MOAB_HAVE_PNETCDF
229 // Wait outside the loop
230 success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
231 if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
232 #endif
233
234 // Read vertices on each local cell, to get localGidVerts and cell connectivity later
235 int verticesOnCellVarId;
236 success = NCFUNC( inq_varid )( _fileId, "elementConn", &verticesOnCellVarId );
237 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of elementConn" );
238 std::vector< int > vertices_on_local_cells( nLocalCells * maxEdgesPerCell );
239 #ifdef MOAB_HAVE_PNETCDF
240 idxReq = 0;
241 #endif
242 indexInArray = 0;
243 for( Range::pair_iterator pair_iter = localGidCells.pair_begin(); pair_iter != localGidCells.pair_end();
244 ++pair_iter )
245 {
246 EntityHandle starth = pair_iter->first;
247 EntityHandle endh = pair_iter->second;
248 NCDF_SIZE read_starts[2] = { static_cast< NCDF_SIZE >( starth - 1 ), 0 };
249 NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( endh - starth + 1 ),
250 static_cast< NCDF_SIZE >( maxEdgesPerCell ) };
251
252 // Do a partial read in each subrange
253 #ifdef MOAB_HAVE_PNETCDF
254 success = NCFUNCREQG( _vara_int )( _fileId, verticesOnCellVarId, read_starts, read_counts,
255 &( vertices_on_local_cells[indexInArray] ), &requests[idxReq++] );
256 #else
257 success = NCFUNCAG( _vara_int )( _fileId, verticesOnCellVarId, read_starts, read_counts,
258 &( vertices_on_local_cells[indexInArray] ) );
259 #endif
260 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read verticesOnCell data in a loop" );
261
262 // Increment the index for next subrange
263 indexInArray += ( endh - starth + 1 ) * maxEdgesPerCell;
264 }
265
266 #ifdef MOAB_HAVE_PNETCDF
267 // Wait outside the loop
268 success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
269 if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
270 #endif
271
272 // Correct local cell vertices array, replace the padded vertices with the last vertices
273 // in the corresponding cells; sometimes the padded vertices are 0, sometimes a large
274 // vertex id. Make sure they are consistent to our padded option
275 for( int local_cell_idx = 0; local_cell_idx < nLocalCells; local_cell_idx++ )
276 {
277 int num_edges = num_edges_on_local_cells[local_cell_idx];
278 int idx_in_local_vert_arr = local_cell_idx * maxEdgesPerCell;
279 int last_vert_idx = vertices_on_local_cells[idx_in_local_vert_arr + num_edges - 1];
280 for( int i = num_edges; i < maxEdgesPerCell; i++ )
281 vertices_on_local_cells[idx_in_local_vert_arr + i] = last_vert_idx;
282 }
283
284 // Create local vertices
285 EntityHandle start_vertex;
286 ErrorCode rval = create_local_vertices( vertices_on_local_cells, start_vertex );MB_CHK_SET_ERR( rval, "Failed to create local vertices for MPAS mesh" );
287
288 // Create local cells, either unpadded or padded
289 if( noMixedElements )
290 {
291 rval = create_padded_local_cells( vertices_on_local_cells, start_vertex, faces );MB_CHK_SET_ERR( rval, "Failed to create padded local cells for MPAS mesh" );
292 }
293 else
294 {
295 rval = create_local_cells( vertices_on_local_cells, num_edges_on_local_cells, start_vertex, faces );MB_CHK_SET_ERR( rval, "Failed to create local cells for MPAS mesh" );
296 }
297
298 return MB_SUCCESS;
299 }
300
301 #ifdef MOAB_HAVE_MPI
302 ErrorCode NCHelperESMF::redistribute_local_cells( int start_cell_idx, ParallelComm* pco )
303 {
304
305 // If possible, apply Zoltan partition
306 // will start from trivial partition in cell space
307 // will read cell connectivities, coordinates of vertices in conn, and then compute centers of the cells
308 //
309 #ifdef MOAB_HAVE_ZOLTAN
310 if( ScdParData::RCBZOLTAN == _readNC->partMethod )
311 {
312 // if we have centerCoords, use them; if not, compute them for Zoltan to work
313 std::vector< double > xverts, yverts, zverts;
314 xverts.resize( nLocalCells );
315 yverts.resize( nLocalCells );
316 zverts.resize( nLocalCells );
317
318 if( centerCoordsId >= 0 )
319 {
320 // read from file
321 std::vector< double > coords( nLocalCells * coordDim );
322
323 NCDF_SIZE read_starts[2] = { static_cast< NCDF_SIZE >( start_cell_idx - 1 ), 0 };
324 NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( nLocalCells ),
325 static_cast< NCDF_SIZE >( coordDim ) };
326 int success = NCFUNCAG( _vara_double )( _fileId, centerCoordsId, read_starts, read_counts, &( coords[0] ) );
327 if( success ) MB_SET_ERR( MB_FAILURE, "Failed on reading center coordinates" );
328 if( 2 == coordDim )
329 {
330 double factor = 1.;
331 if( degrees ) factor = pideg;
332 for( int i = 0; i < nLocalCells; i++ )
333 {
334 double lon = coords[i * 2] * factor;
335 double lat = coords[i * 2 + 1] * factor;
336 double cosphi = cos( lat );
337 zverts[i] = sin( lat );
338 xverts[i] = cosphi * cos( lon );
339 yverts[i] = cosphi * sin( lon );
340 }
341 }
342 }
343 else
344 {
345 // Read connectivities
346 // Read vertices on each local cell, to get localGidVerts and cell connectivity later
347 int verticesOnCellVarId;
348 int success = NCFUNC( inq_varid )( _fileId, "elementConn", &verticesOnCellVarId );
349 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of elementConn" );
350 std::vector< int > vertices_on_local_cells( nLocalCells * maxEdgesPerCell );
351
352 NCDF_SIZE read_starts[2] = { static_cast< NCDF_SIZE >( start_cell_idx - 1 ), 0 };
353 NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( nLocalCells ),
354 static_cast< NCDF_SIZE >( maxEdgesPerCell ) };
355
356 success = NCFUNCAG( _vara_int )( _fileId, verticesOnCellVarId, read_starts, read_counts,
357 &( vertices_on_local_cells[0] ) );
358 if( success ) MB_SET_ERR( MB_FAILURE, "Failed on reading elementConn variable" );
359 std::vector< int > num_edges_on_local_cells( nLocalCells );
360
361 NCDF_SIZE read_start = start_cell_idx - 1;
362 NCDF_SIZE read_count = (NCDF_SIZE)( nLocalCells );
363
364 int nEdgesOnCellVarId;
365 success = NCFUNC( inq_varid )( _fileId, "numElementConn", &nEdgesOnCellVarId );
366 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of numElementConn" );
367 success = NCFUNCAG( _vara_int )( _fileId, nEdgesOnCellVarId, &read_start, &read_count,
368 &( num_edges_on_local_cells[0] ) );
369 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get numElementConn" );
370 // Make a copy of vertices_on_local_cells for sorting (keep original one to set cell
371 // connectivity later)
372
373 // Correct local cell vertices array, replace the padded vertices with the last vertices
374 // in the corresponding cells; sometimes the padded vertices are 0, sometimes a large
375 // vertex id. Make sure they are consistent to our padded option
376 for( int local_cell_idx = 0; local_cell_idx < nLocalCells; local_cell_idx++ )
377 {
378 int num_edges = num_edges_on_local_cells[local_cell_idx];
379 int idx_in_local_vert_arr = local_cell_idx * maxEdgesPerCell;
380 int last_vert_idx = vertices_on_local_cells[idx_in_local_vert_arr + num_edges - 1];
381 for( int i = num_edges; i < maxEdgesPerCell; i++ )
382 vertices_on_local_cells[idx_in_local_vert_arr + i] = last_vert_idx;
383 }
384
385 std::vector< int > vertices_on_local_cells_sorted( vertices_on_local_cells );
386 std::sort( vertices_on_local_cells_sorted.begin(), vertices_on_local_cells_sorted.end() );
387 std::copy( vertices_on_local_cells_sorted.rbegin(), vertices_on_local_cells_sorted.rend(),
388 range_inserter( localGidVerts ) );
389 nLocalVertices = localGidVerts.size();
390
391 // Read nodeCoords for local vertices
392 double* coords = new double[localGidVerts.size() * coordDim];
393 int nodeCoordVarId;
394 success = NCFUNC( inq_varid )( _fileId, "nodeCoords", &nodeCoordVarId );
395 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of nodeCoords" );
396 #ifdef MOAB_HAVE_PNETCDF
397 size_t nb_reads = localGidVerts.psize();
398 std::vector< int > requests( nb_reads );
399 std::vector< int > statuss( nb_reads );
400 size_t idxReq = 0;
401 #endif
402 size_t indexInArray = 0;
403 for( Range::pair_iterator pair_iter = localGidVerts.pair_begin(); pair_iter != localGidVerts.pair_end();
404 ++pair_iter )
405 {
406 EntityHandle starth = pair_iter->first;
407 EntityHandle endh = pair_iter->second;
408 NCDF_SIZE read_starts[2] = { static_cast< NCDF_SIZE >( starth - 1 ), 0 };
409 NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( endh - starth + 1 ),
410 static_cast< NCDF_SIZE >( coordDim ) };
411
412 // Do a partial read in each subrange
413 #ifdef MOAB_HAVE_PNETCDF
414 success = NCFUNCREQG( _vara_double )( _fileId, nodeCoordVarId, read_starts, read_counts,
415 &coords[indexInArray], &requests[idxReq++] );
416 #else
417 success = NCFUNCAG( _vara_double )( _fileId, nodeCoordVarId, read_starts, read_counts,
418 &coords[indexInArray] );
419 #endif
420 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read nodeCoords data in a loop" );
421
422 // Increment the index for next subrange
423 indexInArray += ( endh - starth + 1 ) * coordDim;
424 }
425
426 #ifdef MOAB_HAVE_PNETCDF
427 // Wait outside the loop
428 success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
429 if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
430 #endif
431
432 double* coords3d = new double[localGidVerts.size() * 3];
433 // now convert from lat/lon to 3d
434 if( 2 == coordDim )
435 {
436 // basically convert from degrees to xyz on a sphere
437 double factor = 1.;
438 if( degrees ) factor = pideg;
439 for( int i = 0; i < (int)localGidVerts.size(); i++ )
440 {
441 double lon = coords[i * 2] * factor;
442 double lat = coords[i * 2 + 1] * factor;
443 double cosphi = cos( lat );
444 coords3d[3 * i + 2] = sin( lat );
445 coords3d[3 * i] = cosphi * cos( lon );
446 coords3d[3 * i + 1] = cosphi * sin( lon );
447 }
448 }
449
450 // find the center for each cell
451 for( int i = 0; i < nLocalCells; i++ )
452 {
453 int nv = num_edges_on_local_cells[i];
454 double x = 0., y = 0., z = 0.;
455 for( int j = 0; j < nv; j++ )
456 {
457 int vertexId = vertices_on_local_cells[i * maxEdgesPerCell + j];
458 // now find what index is in gid
459 int index = localGidVerts.index( vertexId ); // it should be a binary search!
460 x += coords3d[3 * index];
461 y += coords3d[3 * index + 1];
462 z += coords3d[3 * index + 2];
463 }
464 if( nv != 0 )
465 {
466 x /= nv;
467 y /= nv;
468 z /= nv;
469 }
470 xverts[i] = x;
471 yverts[i] = y;
472 zverts[i] = z;
473 }
474 delete[] coords3d;
475 }
476 // Zoltan partition using RCB; maybe more studies would be good, as to which partition
477 // is better
478 Interface*& mbImpl = _readNC->mbImpl;
479 DebugOutput& dbgOut = _readNC->dbgOut;
480 ZoltanPartitioner* mbZTool = new ZoltanPartitioner( mbImpl, pco, false, 0, NULL );
481 ErrorCode rval = mbZTool->repartition( xverts, yverts, zverts, start_cell_idx, "RCB", localGidCells );MB_CHK_SET_ERR( rval, "Error in Zoltan partitioning" );
482 delete mbZTool;
483
484 dbgOut.tprintf( 1, "After Zoltan partitioning, localGidCells.psize() = %d\n", (int)localGidCells.psize() );
485 dbgOut.tprintf( 1, " localGidCells.size() = %d\n", (int)localGidCells.size() );
486
487 // This is important: local cells are now redistributed, so nLocalCells might be different!
488 nLocalCells = localGidCells.size();
489 localGidVerts.clear();
490 return MB_SUCCESS;
491 }
492 #endif
493
494 // By default, apply trivial partition
495 localGidCells.insert( start_cell_idx, start_cell_idx + nLocalCells - 1 );
496
497 return MB_SUCCESS;
498 }
499 #endif
500
501 ErrorCode NCHelperESMF::create_local_vertices( const std::vector< int >& vertices_on_local_cells,
502 EntityHandle& start_vertex )
503 {
504 Interface*& mbImpl = _readNC->mbImpl;
505 Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
506 const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
507 DebugOutput& dbgOut = _readNC->dbgOut;
508
509 // Make a copy of vertices_on_local_cells for sorting (keep original one to set cell
510 // connectivity later)
511 std::vector< int > vertices_on_local_cells_sorted( vertices_on_local_cells );
512 std::sort( vertices_on_local_cells_sorted.begin(), vertices_on_local_cells_sorted.end() );
513 std::copy( vertices_on_local_cells_sorted.rbegin(), vertices_on_local_cells_sorted.rend(),
514 range_inserter( localGidVerts ) );
515 nLocalVertices = localGidVerts.size();
516
517 dbgOut.tprintf( 1, " localGidVerts.psize() = %d\n", (int)localGidVerts.psize() );
518 dbgOut.tprintf( 1, " localGidVerts.size() = %d\n", (int)localGidVerts.size() );
519
520 // Create local vertices
521 std::vector< double* > arrays;
522 ErrorCode rval =
523 _readNC->readMeshIface->get_node_coords( 3, nLocalVertices, 0, start_vertex, arrays, nLocalVertices );MB_CHK_SET_ERR( rval, "Failed to create local vertices" );
524
525 // Add local vertices to current file set
526 Range local_verts_range( start_vertex, start_vertex + nLocalVertices - 1 );
527 rval = _readNC->mbImpl->add_entities( _fileSet, local_verts_range );MB_CHK_SET_ERR( rval, "Failed to add local vertices to current file set" );
528
529 // Get ptr to GID memory for local vertices
530 int count = 0;
531 void* data = NULL;
532 rval = mbImpl->tag_iterate( mGlobalIdTag, local_verts_range.begin(), local_verts_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate global id tag on local vertices" );
533 assert( count == nLocalVertices );
534 int* gid_data = (int*)data;
535 std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
536
537 // Duplicate GID data, which will be used to resolve sharing
538 if( mpFileIdTag )
539 {
540 rval = mbImpl->tag_iterate( *mpFileIdTag, local_verts_range.begin(), local_verts_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate file id tag on local vertices" );
541 assert( count == nLocalVertices );
542 int bytes_per_tag = 4;
543 rval = mbImpl->tag_get_bytes( *mpFileIdTag, bytes_per_tag );MB_CHK_SET_ERR( rval, "Can't get number of bytes for file id tag" );
544 if( 4 == bytes_per_tag )
545 {
546 gid_data = (int*)data;
547 std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
548 }
549 else if( 8 == bytes_per_tag )
550 { // Should be a handle tag on 64 bit machine?
551 long* handle_tag_data = (long*)data;
552 std::copy( localGidVerts.begin(), localGidVerts.end(), handle_tag_data );
553 }
554 }
555
556 #ifdef MOAB_HAVE_PNETCDF
557 size_t nb_reads = localGidVerts.psize();
558 std::vector< int > requests( nb_reads );
559 std::vector< int > statuss( nb_reads );
560 size_t idxReq = 0;
561 #endif
562
563 // Read nodeCoords for local vertices
564 double* coords = new double[localGidVerts.size() * coordDim];
565 int nodeCoordVarId;
566 int success = NCFUNC( inq_varid )( _fileId, "nodeCoords", &nodeCoordVarId );
567 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of nodeCoords" );
568 size_t indexInArray = 0;
569 for( Range::pair_iterator pair_iter = localGidVerts.pair_begin(); pair_iter != localGidVerts.pair_end();
570 ++pair_iter )
571 {
572 EntityHandle starth = pair_iter->first;
573 EntityHandle endh = pair_iter->second;
574 NCDF_SIZE read_starts[2] = { static_cast< NCDF_SIZE >( starth - 1 ), 0 };
575 NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( endh - starth + 1 ),
576 static_cast< NCDF_SIZE >( coordDim ) };
577
578 // Do a partial read in each subrange
579 #ifdef MOAB_HAVE_PNETCDF
580 success = NCFUNCREQG( _vara_double )( _fileId, nodeCoordVarId, read_starts, read_counts, &coords[indexInArray],
581 &requests[idxReq++] );
582 #else
583 success = NCFUNCAG( _vara_double )( _fileId, nodeCoordVarId, read_starts, read_counts, &coords[indexInArray] );
584 #endif
585 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read nodeCoords data in a loop" );
586
587 // Increment the index for next subrange
588 indexInArray += ( endh - starth + 1 ) * coordDim;
589 }
590
591 #ifdef MOAB_HAVE_PNETCDF
592 // Wait outside the loop
593 success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
594 if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
595 #endif
596
597 // now convert from lat/lon to 3d
598 if( 2 == coordDim )
599 {
600 // basically convert from degrees to xyz on a sphere
601 double factor = 1;
602 if( degrees ) factor = pideg;
603 for( int i = 0; i < (int)localGidVerts.size(); i++ )
604 {
605 double lon = coords[i * 2] * factor;
606 double lat = coords[i * 2 + 1] * factor;
607 double cosphi = cos( lat );
608 arrays[2][i] = sin( lat );
609 arrays[0][i] = cosphi * cos( lon );
610 arrays[1][i] = cosphi * sin( lon );
611 }
612 }
613 delete[] coords;
614 return MB_SUCCESS;
615 }
616
617 ErrorCode NCHelperESMF::create_local_cells( const std::vector< int >& vertices_on_local_cells,
618 const std::vector< int >& num_edges_on_local_cells,
619 EntityHandle start_vertex,
620 Range& faces )
621 {
622 Interface*& mbImpl = _readNC->mbImpl;
623 Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
624
625 // Divide local cells into groups based on the number of edges
626 Range local_cells_with_n_edges[DEFAULT_MAX_EDGES_PER_CELL + 1];
627 // Insert larger values before smaller ones to increase efficiency
628 for( int i = nLocalCells - 1; i >= 0; i-- )
629 {
630 int num_edges = num_edges_on_local_cells[i];
631 local_cells_with_n_edges[num_edges].insert( localGidCells[i] ); // Global cell index
632 }
633
634 std::vector< int > num_edges_on_cell_groups;
635 for( int i = 3; i <= maxEdgesPerCell; i++ )
636 {
637 if( local_cells_with_n_edges[i].size() > 0 ) num_edges_on_cell_groups.push_back( i );
638 }
639 int numCellGroups = (int)num_edges_on_cell_groups.size();
640 EntityHandle* conn_arr_local_cells_with_n_edges[DEFAULT_MAX_EDGES_PER_CELL + 1];
641 for( int i = 0; i < numCellGroups; i++ )
642 {
643 int num_edges_per_cell = num_edges_on_cell_groups[i];
644 int num_group_cells = (int)local_cells_with_n_edges[num_edges_per_cell].size();
645
646 EntityType typeEl = MBTRI;
647 if( num_edges_per_cell == 4 ) typeEl = MBQUAD;
648 if( num_edges_per_cell > 4 ) typeEl = MBPOLYGON;
649 // Create local cells for each non-empty cell group
650 EntityHandle start_element;
651 ErrorCode rval =
652 _readNC->readMeshIface->get_element_connect( num_group_cells, num_edges_per_cell, typeEl, 0, start_element,
653 conn_arr_local_cells_with_n_edges[num_edges_per_cell],
654 num_group_cells );MB_CHK_SET_ERR( rval, "Failed to create local cells" );
655 faces.insert( start_element, start_element + num_group_cells - 1 );
656
657 // Add local cells to current file set
658 Range local_cells_range( start_element, start_element + num_group_cells - 1 );
659 rval = _readNC->mbImpl->add_entities( _fileSet, local_cells_range );MB_CHK_SET_ERR( rval, "Failed to add local cells to current file set" );
660
661 // Get ptr to gid memory for local cells
662 int count = 0;
663 void* data = NULL;
664 rval = mbImpl->tag_iterate( mGlobalIdTag, local_cells_range.begin(), local_cells_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate global id tag on local cells" );
665 assert( count == num_group_cells );
666 int* gid_data = (int*)data;
667 std::copy( local_cells_with_n_edges[num_edges_per_cell].begin(),
668 local_cells_with_n_edges[num_edges_per_cell].end(), gid_data );
669
670 // Set connectivity array with proper local vertices handles
671 for( int j = 0; j < num_group_cells; j++ )
672 {
673 EntityHandle global_cell_idx =
674 local_cells_with_n_edges[num_edges_per_cell][j]; // Global cell index, 1 based
675 int local_cell_idx = localGidCells.index( global_cell_idx ); // Local cell index, 0 based
676 assert( local_cell_idx != -1 );
677
678 for( int k = 0; k < num_edges_per_cell; k++ )
679 {
680 EntityHandle global_vert_idx =
681 vertices_on_local_cells[local_cell_idx * maxEdgesPerCell + k]; // Global vertex index, 1 based
682 int local_vert_idx = localGidVerts.index( global_vert_idx ); // Local vertex index, 0 based
683 assert( local_vert_idx != -1 );
684 conn_arr_local_cells_with_n_edges[num_edges_per_cell][j * num_edges_per_cell + k] =
685 start_vertex + local_vert_idx;
686 }
687 }
688 }
689
690 return MB_SUCCESS;
691 }
692
693 ErrorCode NCHelperESMF::create_padded_local_cells( const std::vector< int >& vertices_on_local_cells,
694 EntityHandle start_vertex,
695 Range& faces )
696 {
697 Interface*& mbImpl = _readNC->mbImpl;
698 Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
699
700 // Create cells for this cell group
701 EntityHandle start_element;
702 EntityHandle* conn_arr_local_cells = NULL;
703 ErrorCode rval = _readNC->readMeshIface->get_element_connect( nLocalCells, maxEdgesPerCell, MBPOLYGON, 0,
704 start_element, conn_arr_local_cells, nLocalCells );MB_CHK_SET_ERR( rval, "Failed to create local cells" );
705 faces.insert( start_element, start_element + nLocalCells - 1 );
706
707 // Add local cells to current file set
708 Range local_cells_range( start_element, start_element + nLocalCells - 1 );
709 rval = _readNC->mbImpl->add_entities( _fileSet, local_cells_range );MB_CHK_SET_ERR( rval, "Failed to add local cells to current file set" );
710
711 // Get ptr to GID memory for local cells
712 int count = 0;
713 void* data = NULL;
714 rval = mbImpl->tag_iterate( mGlobalIdTag, local_cells_range.begin(), local_cells_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate global id tag on local cells" );
715 assert( count == nLocalCells );
716 int* gid_data = (int*)data;
717 std::copy( localGidCells.begin(), localGidCells.end(), gid_data );
718
719 // Set connectivity array with proper local vertices handles
720 // vertices_on_local_cells array was already corrected to have the last vertices padded
721 // no need for extra checks considering
722 for( int local_cell_idx = 0; local_cell_idx < nLocalCells; local_cell_idx++ )
723 {
724 for( int i = 0; i < maxEdgesPerCell; i++ )
725 {
726 EntityHandle global_vert_idx =
727 vertices_on_local_cells[local_cell_idx * maxEdgesPerCell + i]; // Global vertex index, 1 based
728 int local_vert_idx = localGidVerts.index( global_vert_idx ); // Local vertex index, 0 based
729 assert( local_vert_idx != -1 );
730 conn_arr_local_cells[local_cell_idx * maxEdgesPerCell + i] = start_vertex + local_vert_idx;
731 }
732 }
733
734 return MB_SUCCESS;
735 }
736
737 } // namespace moab
1.9.1.