Branch data Line data Source code
1 : : #include "NCHelperGCRM.hpp"
2 : : #include "moab/ReadUtilIface.hpp"
3 : : #include "moab/FileOptions.hpp"
4 : : #include "moab/SpectralMeshTool.hpp"
5 : : #include "MBTagConventions.hpp"
6 : :
7 : : #ifdef MOAB_HAVE_ZOLTAN
8 : : #include "moab/ZoltanPartitioner.hpp"
9 : : #endif
10 : :
11 : : #include <cmath>
12 : :
13 : : namespace moab
14 : : {
15 : :
16 : : // GCRM cells are either pentagons or hexagons, and pentagons are always padded to hexagons
17 : : const int EDGES_PER_CELL = 6;
18 : :
19 : 0 : NCHelperGCRM::NCHelperGCRM( ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet )
20 [ # # ]: 0 : : UcdNCHelper( readNC, fileId, opts, fileSet ), createGatherSet( false )
21 : : {
22 : : // Ignore variables containing topological information
23 [ # # ][ # # ]: 0 : ignoredVarNames.insert( "grid" );
24 [ # # ][ # # ]: 0 : ignoredVarNames.insert( "cell_corners" );
25 [ # # ][ # # ]: 0 : ignoredVarNames.insert( "cell_edges" );
26 [ # # ][ # # ]: 0 : ignoredVarNames.insert( "edge_corners" );
27 [ # # ][ # # ]: 0 : ignoredVarNames.insert( "cell_neighbors" );
28 : 0 : }
29 : :
30 : 0 : bool NCHelperGCRM::can_read_file( ReadNC* readNC )
31 : : {
32 : 0 : std::vector< std::string >& dimNames = readNC->dimNames;
33 : :
34 : : // If dimension name "cells" exists then it should be the GCRM grid
35 [ # # ][ # # ]: 0 : if( std::find( dimNames.begin(), dimNames.end(), std::string( "cells" ) ) != dimNames.end() ) return true;
[ # # ][ # # ]
36 : :
37 : 0 : return false;
38 : : }
39 : :
40 : 0 : ErrorCode NCHelperGCRM::init_mesh_vals()
41 : : {
42 : 0 : std::vector< std::string >& dimNames = _readNC->dimNames;
43 : 0 : std::vector< int >& dimLens = _readNC->dimLens;
44 : 0 : std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
45 : :
46 : : ErrorCode rval;
47 : : unsigned int idx;
48 : 0 : std::vector< std::string >::iterator vit;
49 : :
50 : : // Look for time dimension
51 [ # # ][ # # ]: 0 : if( ( vit = std::find( dimNames.begin(), dimNames.end(), "Time" ) ) != dimNames.end() )
[ # # ]
52 [ # # ]: 0 : idx = vit - dimNames.begin();
53 [ # # ][ # # ]: 0 : else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
[ # # ]
54 [ # # ]: 0 : idx = vit - dimNames.begin();
55 : : else
56 : : {
57 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't find 'Time' or 'time' dimension" );
[ # # ][ # # ]
[ # # ]
58 : : }
59 : 0 : tDim = idx;
60 [ # # ]: 0 : nTimeSteps = dimLens[idx];
61 : :
62 : : // Get number of cells
63 [ # # ][ # # ]: 0 : if( ( vit = std::find( dimNames.begin(), dimNames.end(), "cells" ) ) != dimNames.end() )
[ # # ]
64 [ # # ]: 0 : idx = vit - dimNames.begin();
65 : : else
66 : : {
67 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't find 'cells' dimension" );
[ # # ][ # # ]
[ # # ]
68 : : }
69 : 0 : cDim = idx;
70 [ # # ]: 0 : nCells = dimLens[idx];
71 : :
72 : : // Get number of edges
73 [ # # ][ # # ]: 0 : if( ( vit = std::find( dimNames.begin(), dimNames.end(), "edges" ) ) != dimNames.end() )
[ # # ]
74 [ # # ]: 0 : idx = vit - dimNames.begin();
75 : : else
76 : : {
77 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't find 'edges' dimension" );
[ # # ][ # # ]
[ # # ]
78 : : }
79 : 0 : eDim = idx;
80 [ # # ]: 0 : nEdges = dimLens[idx];
81 : :
82 : : // Get number of vertices
83 : 0 : vDim = -1;
84 [ # # ][ # # ]: 0 : if( ( vit = std::find( dimNames.begin(), dimNames.end(), "corners" ) ) != dimNames.end() )
[ # # ]
85 [ # # ]: 0 : idx = vit - dimNames.begin();
86 : : else
87 : : {
88 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't find 'corners' dimension" );
[ # # ][ # # ]
[ # # ]
89 : : }
90 : 0 : vDim = idx;
91 [ # # ]: 0 : nVertices = dimLens[idx];
92 : :
93 : : // Get number of layers
94 [ # # ][ # # ]: 0 : if( ( vit = std::find( dimNames.begin(), dimNames.end(), "layers" ) ) != dimNames.end() )
[ # # ]
95 [ # # ]: 0 : idx = vit - dimNames.begin();
96 : : else
97 : : {
98 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't find 'layers' dimension" );
[ # # ][ # # ]
[ # # ]
99 : : }
100 : 0 : levDim = idx;
101 [ # # ]: 0 : nLevels = dimLens[idx];
102 : :
103 : : // Dimension indices for other optional levels
104 [ # # ]: 0 : std::vector< unsigned int > opt_lev_dims;
105 : :
106 : : // Get index of interface levels
107 [ # # ][ # # ]: 0 : if( ( vit = std::find( dimNames.begin(), dimNames.end(), "interfaces" ) ) != dimNames.end() )
[ # # ]
108 : : {
109 [ # # ]: 0 : idx = vit - dimNames.begin();
110 [ # # ]: 0 : opt_lev_dims.push_back( idx );
111 : : }
112 : :
113 [ # # ]: 0 : std::map< std::string, ReadNC::VarData >::iterator vmit;
114 : :
115 : : // Store time coordinate values in tVals
116 [ # # ]: 0 : if( nTimeSteps > 0 )
117 : : {
118 [ # # ][ # # ]: 0 : if( ( vmit = varInfo.find( "time" ) ) != varInfo.end() )
[ # # ][ # # ]
119 : : {
120 [ # # ][ # # ]: 0 : rval = read_coordinate( "time", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 'time' variable" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
121 : : }
122 [ # # ][ # # ]: 0 : else if( ( vmit = varInfo.find( "t" ) ) != varInfo.end() )
[ # # ][ # # ]
123 : : {
124 [ # # ][ # # ]: 0 : rval = read_coordinate( "t", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 't' variable" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
125 : : }
126 : : else
127 : : {
128 : : // If expected time variable is not available, set dummy time coordinate values to tVals
129 [ # # ]: 0 : for( int t = 0; t < nTimeSteps; t++ )
130 [ # # ]: 0 : tVals.push_back( (double)t );
131 : : }
132 : : }
133 : :
134 : : // For each variable, determine the entity location type and number of levels
135 [ # # ][ # # ]: 0 : for( vmit = varInfo.begin(); vmit != varInfo.end(); ++vmit )
[ # # ]
136 : : {
137 [ # # ]: 0 : ReadNC::VarData& vd = ( *vmit ).second;
138 : :
139 : : // Default entLoc is ENTLOCSET
140 [ # # ][ # # ]: 0 : if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
[ # # ]
141 : : {
142 [ # # ][ # # ]: 0 : if( std::find( vd.varDims.begin(), vd.varDims.end(), vDim ) != vd.varDims.end() )
[ # # ]
143 : 0 : vd.entLoc = ReadNC::ENTLOCVERT;
144 [ # # ][ # # ]: 0 : else if( std::find( vd.varDims.begin(), vd.varDims.end(), eDim ) != vd.varDims.end() )
[ # # ]
145 : 0 : vd.entLoc = ReadNC::ENTLOCEDGE;
146 [ # # ][ # # ]: 0 : else if( std::find( vd.varDims.begin(), vd.varDims.end(), cDim ) != vd.varDims.end() )
[ # # ]
147 : 0 : vd.entLoc = ReadNC::ENTLOCFACE;
148 : : }
149 : :
150 : : // Default numLev is 0
151 [ # # ][ # # ]: 0 : if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() )
[ # # ]
152 : 0 : vd.numLev = nLevels;
153 : : else
154 : : {
155 : : // If layers dimension is not found, try other optional levels such as interfaces
156 [ # # ]: 0 : for( unsigned int i = 0; i < opt_lev_dims.size(); i++ )
157 : : {
158 [ # # ][ # # ]: 0 : if( std::find( vd.varDims.begin(), vd.varDims.end(), opt_lev_dims[i] ) != vd.varDims.end() )
[ # # ][ # # ]
159 : : {
160 [ # # ][ # # ]: 0 : vd.numLev = dimLens[opt_lev_dims[i]];
161 : 0 : break;
162 : : }
163 : : }
164 : : }
165 : : }
166 : :
167 : : // Hack: create dummy variables for dimensions (like cells) with no corresponding coordinate
168 : : // variables
169 [ # # ][ # # ]: 0 : rval = create_dummy_variables();MB_CHK_SET_ERR( rval, "Failed to create dummy variables" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
170 : :
171 : 0 : return MB_SUCCESS;
172 : : }
173 : :
174 : : // When noMesh option is used on this read, the old ReadNC class instance for last read can get out
175 : : // of scope (and deleted). The old instance initialized some variables properly when the mesh was
176 : : // created, but they are now lost. The new instance (will not create the mesh with noMesh option)
177 : : // has to restore them based on the existing mesh from last read
178 : 0 : ErrorCode NCHelperGCRM::check_existing_mesh()
179 : : {
180 : 0 : Interface*& mbImpl = _readNC->mbImpl;
181 : 0 : Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
182 : 0 : bool& noMesh = _readNC->noMesh;
183 : :
184 [ # # ]: 0 : if( noMesh )
185 : : {
186 : : ErrorCode rval;
187 : :
188 [ # # ]: 0 : if( localGidVerts.empty() )
189 : : {
190 : : // Get all vertices from current file set (it is the input set in no_mesh scenario)
191 [ # # ]: 0 : Range local_verts;
192 [ # # ][ # # ]: 0 : rval = mbImpl->get_entities_by_dimension( _fileSet, 0, local_verts );MB_CHK_SET_ERR( rval, "Trouble getting local vertices in current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
193 : :
194 [ # # ][ # # ]: 0 : if( !local_verts.empty() )
195 : : {
196 [ # # ][ # # ]: 0 : std::vector< int > gids( local_verts.size() );
197 : :
198 : : // !IMPORTANT : this has to be the GLOBAL_ID tag
199 [ # # ][ # # ]: 0 : rval = mbImpl->tag_get_data( mGlobalIdTag, local_verts, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting local gid values of vertices" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
200 : :
201 : : // Restore localGidVerts
202 [ # # ][ # # ]: 0 : std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidVerts ) );
203 [ # # ][ # # ]: 0 : nLocalVertices = localGidVerts.size();
[ # # ]
204 : 0 : }
205 : : }
206 : :
207 [ # # ]: 0 : if( localGidEdges.empty() )
208 : : {
209 : : // Get all edges from current file set (it is the input set in no_mesh scenario)
210 [ # # ]: 0 : Range local_edges;
211 [ # # ][ # # ]: 0 : rval = mbImpl->get_entities_by_dimension( _fileSet, 1, local_edges );MB_CHK_SET_ERR( rval, "Trouble getting local edges in current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
212 : :
213 [ # # ][ # # ]: 0 : if( !local_edges.empty() )
214 : : {
215 [ # # ][ # # ]: 0 : std::vector< int > gids( local_edges.size() );
216 : :
217 : : // !IMPORTANT : this has to be the GLOBAL_ID tag
218 [ # # ][ # # ]: 0 : rval = mbImpl->tag_get_data( mGlobalIdTag, local_edges, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting local gid values of edges" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
219 : :
220 : : // Restore localGidEdges
221 [ # # ][ # # ]: 0 : std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidEdges ) );
222 [ # # ][ # # ]: 0 : nLocalEdges = localGidEdges.size();
[ # # ]
223 : 0 : }
224 : : }
225 : :
226 [ # # ]: 0 : if( localGidCells.empty() )
227 : : {
228 : : // Get all cells from current file set (it is the input set in no_mesh scenario)
229 [ # # ]: 0 : Range local_cells;
230 [ # # ][ # # ]: 0 : rval = mbImpl->get_entities_by_dimension( _fileSet, 2, local_cells );MB_CHK_SET_ERR( rval, "Trouble getting local cells in current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
231 : :
232 [ # # ][ # # ]: 0 : if( !local_cells.empty() )
233 : : {
234 [ # # ][ # # ]: 0 : std::vector< int > gids( local_cells.size() );
235 : :
236 : : // !IMPORTANT : this has to be the GLOBAL_ID tag
237 [ # # ][ # # ]: 0 : rval = mbImpl->tag_get_data( mGlobalIdTag, local_cells, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting local gid values of cells" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
238 : :
239 : : // Restore localGidCells
240 [ # # ][ # # ]: 0 : std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidCells ) );
241 [ # # ][ # # ]: 0 : nLocalCells = localGidCells.size();
[ # # ]
242 : 0 : }
243 : : }
244 : : }
245 : :
246 : 0 : return MB_SUCCESS;
247 : : }
248 : :
249 : 0 : ErrorCode NCHelperGCRM::create_mesh( Range& faces )
250 : : {
251 : 0 : bool& noEdges = _readNC->noEdges;
252 : 0 : DebugOutput& dbgOut = _readNC->dbgOut;
253 : :
254 : : #ifdef MOAB_HAVE_MPI
255 : 0 : int rank = 0;
256 : 0 : int procs = 1;
257 : 0 : bool& isParallel = _readNC->isParallel;
258 [ # # ]: 0 : if( isParallel )
259 : : {
260 : 0 : ParallelComm*& myPcomm = _readNC->myPcomm;
261 [ # # ][ # # ]: 0 : rank = myPcomm->proc_config().proc_rank();
262 [ # # ][ # # ]: 0 : procs = myPcomm->proc_config().proc_size();
263 : : }
264 : :
265 : : // Need to know whether we'll be creating gather mesh
266 [ # # ]: 0 : if( rank == _readNC->gatherSetRank ) createGatherSet = true;
267 : :
268 [ # # ]: 0 : if( procs >= 2 )
269 : : {
270 : : // Shift rank to obtain a rotated trivial partition
271 : 0 : int shifted_rank = rank;
272 : 0 : int& trivialPartitionShift = _readNC->trivialPartitionShift;
273 [ # # ]: 0 : if( trivialPartitionShift > 0 ) shifted_rank = ( rank + trivialPartitionShift ) % procs;
274 : :
275 : : // Compute the number of local cells on this proc
276 : 0 : nLocalCells = int( std::floor( 1.0 * nCells / procs ) );
277 : :
278 : : // The starting global cell index in the GCRM file for this proc
279 : 0 : int start_cell_idx = shifted_rank * nLocalCells;
280 : :
281 : : // Number of extra cells after equal split over procs
282 : 0 : int iextra = nCells % procs;
283 : :
284 : : // Allocate extra cells over procs
285 [ # # ]: 0 : if( shifted_rank < iextra ) nLocalCells++;
286 [ # # ]: 0 : start_cell_idx += std::min( shifted_rank, iextra );
287 : :
288 : 0 : start_cell_idx++; // 0 based -> 1 based
289 : :
290 : : // Redistribute local cells after trivial partition (e.g. apply Zoltan partition)
291 [ # # ][ # # ]: 0 : ErrorCode rval = redistribute_local_cells( start_cell_idx );MB_CHK_SET_ERR( rval, "Failed to redistribute local cells after trivial partition" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
292 : : }
293 : : else
294 : : {
295 : 0 : nLocalCells = nCells;
296 [ # # ]: 0 : localGidCells.insert( 1, nLocalCells );
297 : : }
298 : : #else
299 : : nLocalCells = nCells;
300 : : localGidCells.insert( 1, nLocalCells );
301 : : #endif
302 [ # # ][ # # ]: 0 : dbgOut.tprintf( 1, " localGidCells.psize() = %d\n", (int)localGidCells.psize() );
303 [ # # ][ # # ]: 0 : dbgOut.tprintf( 1, " localGidCells.size() = %d\n", (int)localGidCells.size() );
304 : :
305 : : // Read vertices on each local cell, to get localGidVerts and cell connectivity later
306 : : int verticesOnCellVarId;
307 [ # # ]: 0 : int success = NCFUNC( inq_varid )( _fileId, "cell_corners", &verticesOnCellVarId );
308 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of cell_corners" );
[ # # ][ # # ]
[ # # ][ # # ]
309 [ # # ]: 0 : std::vector< int > vertices_on_local_cells( nLocalCells * EDGES_PER_CELL );
310 [ # # ]: 0 : dbgOut.tprintf( 1, " nLocalCells = %d\n", (int)nLocalCells );
311 [ # # ]: 0 : dbgOut.tprintf( 1, " vertices_on_local_cells.size() = %d\n", (int)vertices_on_local_cells.size() );
312 : : #ifdef MOAB_HAVE_PNETCDF
313 : : size_t nb_reads = localGidCells.psize();
314 : : std::vector< int > requests( nb_reads );
315 : : std::vector< int > statuss( nb_reads );
316 : : size_t idxReq = 0;
317 : : #endif
318 : 0 : size_t indexInArray = 0;
319 [ # # ][ # # ]: 0 : for( Range::pair_iterator pair_iter = localGidCells.pair_begin(); pair_iter != localGidCells.pair_end();
[ # # ][ # # ]
[ # # ]
320 : : ++pair_iter )
321 : : {
322 [ # # ]: 0 : EntityHandle starth = pair_iter->first;
323 [ # # ]: 0 : EntityHandle endh = pair_iter->second;
324 [ # # ]: 0 : dbgOut.tprintf( 1, " cell_corners starth = %d\n", (int)starth );
325 [ # # ]: 0 : dbgOut.tprintf( 1, " cell_corners endh = %d\n", (int)endh );
326 : 0 : NCDF_SIZE read_starts[2] = { static_cast< NCDF_SIZE >( starth - 1 ), 0 };
327 : 0 : NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( endh - starth + 1 ),
328 : 0 : static_cast< NCDF_SIZE >( EDGES_PER_CELL ) };
329 : :
330 : : // Do a partial read in each subrange
331 : : #ifdef MOAB_HAVE_PNETCDF
332 : : success = NCFUNCREQG( _vara_int )( _fileId, verticesOnCellVarId, read_starts, read_counts,
333 : : &( vertices_on_local_cells[indexInArray] ), &requests[idxReq++] );
334 : : #else
335 : : success = NCFUNCAG( _vara_int )( _fileId, verticesOnCellVarId, read_starts, read_counts,
336 [ # # ][ # # ]: 0 : &( vertices_on_local_cells[indexInArray] ) );
337 : : #endif
338 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read cell_corners data in a loop" );
[ # # ][ # # ]
[ # # ][ # # ]
339 : :
340 : : // Increment the index for next subrange
341 : 0 : indexInArray += ( endh - starth + 1 ) * EDGES_PER_CELL;
342 : : }
343 : :
344 : : #ifdef MOAB_HAVE_PNETCDF
345 : : // Wait outside the loop
346 : : success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
347 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
348 : : #endif
349 : :
350 : : // Correct vertices_on_local_cells array. Pentagons as hexagons should have
351 : : // a connectivity like 123455 and not 122345
352 [ # # ]: 0 : for( int local_cell_idx = 0; local_cell_idx < nLocalCells; local_cell_idx++ )
353 : : {
354 [ # # ]: 0 : int* pvertex = &vertices_on_local_cells[local_cell_idx * EDGES_PER_CELL];
355 [ # # ]: 0 : for( int k = 0; k < EDGES_PER_CELL - 2; k++ )
356 : : {
357 [ # # ]: 0 : if( *( pvertex + k ) == *( pvertex + k + 1 ) )
358 : : {
359 : : // Shift the connectivity
360 [ # # ]: 0 : for( int kk = k + 1; kk < EDGES_PER_CELL - 1; kk++ )
361 : 0 : *( pvertex + kk ) = *( pvertex + kk + 1 );
362 : : // No need to try next k
363 : 0 : break;
364 : : }
365 : : }
366 : : }
367 : :
368 : : // GCRM is 0 based, convert vertex indices from 0 to 1 based
369 [ # # ]: 0 : for( std::size_t idx = 0; idx < vertices_on_local_cells.size(); idx++ )
370 [ # # ]: 0 : vertices_on_local_cells[idx] += 1;
371 : :
372 : : // Create local vertices
373 : : EntityHandle start_vertex;
374 [ # # ][ # # ]: 0 : ErrorCode rval = create_local_vertices( vertices_on_local_cells, start_vertex );MB_CHK_SET_ERR( rval, "Failed to create local vertices for GCRM mesh" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
375 : :
376 : : // Create local edges (unless NO_EDGES read option is set)
377 [ # # ]: 0 : if( !noEdges )
378 : : {
379 [ # # ][ # # ]: 0 : rval = create_local_edges( start_vertex );MB_CHK_SET_ERR( rval, "Failed to create local edges for GCRM mesh" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
380 : : }
381 : :
382 : : // Create local cells, padded
383 [ # # ][ # # ]: 0 : rval = create_padded_local_cells( vertices_on_local_cells, start_vertex, faces );MB_CHK_SET_ERR( rval, "Failed to create padded local cells for GCRM mesh" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
384 : :
385 [ # # ]: 0 : if( createGatherSet )
386 : : {
387 : : EntityHandle gather_set;
388 [ # # ][ # # ]: 0 : rval = _readNC->readMeshIface->create_gather_set( gather_set );MB_CHK_SET_ERR( rval, "Failed to create gather set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
389 : :
390 : : // Create gather set vertices
391 : : EntityHandle start_gather_set_vertex;
392 [ # # ][ # # ]: 0 : rval = create_gather_set_vertices( gather_set, start_gather_set_vertex );MB_CHK_SET_ERR( rval, "Failed to create gather set vertices for GCRM mesh" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
393 : :
394 : : // Create gather set edges (unless NO_EDGES read option is set)
395 [ # # ]: 0 : if( !noEdges )
396 : : {
397 [ # # ][ # # ]: 0 : rval = create_gather_set_edges( gather_set, start_gather_set_vertex );MB_CHK_SET_ERR( rval, "Failed to create gather set edges for GCRM mesh" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
398 : : }
399 : :
400 : : // Create gather set cells with padding
401 [ # # ][ # # ]: 0 : rval = create_padded_gather_set_cells( gather_set, start_gather_set_vertex );MB_CHK_SET_ERR( rval, "Failed to create padded gather set cells for GCRM mesh" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
402 : : }
403 : :
404 : 0 : return MB_SUCCESS;
405 : : }
406 : :
407 : 0 : ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset_allocate( std::vector< ReadNC::VarData >& vdatas,
408 : : std::vector< int >& tstep_nums )
409 : : {
410 : 0 : Interface*& mbImpl = _readNC->mbImpl;
411 : 0 : std::vector< int >& dimLens = _readNC->dimLens;
412 : 0 : bool& noEdges = _readNC->noEdges;
413 : 0 : DebugOutput& dbgOut = _readNC->dbgOut;
414 : :
415 : 0 : Range* range = NULL;
416 : :
417 : : // Get vertices
418 [ # # ]: 0 : Range verts;
419 [ # # ][ # # ]: 0 : ErrorCode rval = mbImpl->get_entities_by_dimension( _fileSet, 0, verts );MB_CHK_SET_ERR( rval, "Trouble getting vertices in current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
420 [ # # ][ # # ]: 0 : assert( "Should only have a single vertex subrange, since they were read in one shot" && verts.psize() == 1 );
421 : :
422 : : // Get edges
423 [ # # ]: 0 : Range edges;
424 [ # # ][ # # ]: 0 : rval = mbImpl->get_entities_by_dimension( _fileSet, 1, edges );MB_CHK_SET_ERR( rval, "Trouble getting edges in current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
425 : :
426 : : // Get faces
427 [ # # ]: 0 : Range faces;
428 [ # # ][ # # ]: 0 : rval = mbImpl->get_entities_by_dimension( _fileSet, 2, faces );MB_CHK_SET_ERR( rval, "Trouble getting faces in current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
429 [ # # ][ # # ]: 0 : assert( "Should only have a single face subrange, since they were read in one shot" && faces.psize() == 1 );
430 : :
431 : : #ifdef MOAB_HAVE_MPI
432 : 0 : bool& isParallel = _readNC->isParallel;
433 [ # # ]: 0 : if( isParallel )
434 : : {
435 : 0 : ParallelComm*& myPcomm = _readNC->myPcomm;
436 [ # # ][ # # ]: 0 : rval = myPcomm->filter_pstatus( faces, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &facesOwned );MB_CHK_SET_ERR( rval, "Trouble getting owned faces in current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
437 : : }
438 : : else
439 [ # # ]: 0 : facesOwned = faces; // not running in parallel, but still with MPI
440 : : #else
441 : : facesOwned = faces;
442 : : #endif
443 : :
444 [ # # ]: 0 : for( unsigned int i = 0; i < vdatas.size(); i++ )
445 : : {
446 : : // Skip edge variables, if specified by the read options
447 [ # # ][ # # ]: 0 : if( noEdges && vdatas[i].entLoc == ReadNC::ENTLOCEDGE ) continue;
[ # # ][ # # ]
448 : :
449 : : // Support non-set variables with 3 dimensions like (time, cells, layers), or
450 : : // 2 dimensions like (time, cells)
451 [ # # ][ # # ]: 0 : assert( 3 == vdatas[i].varDims.size() || 2 == vdatas[i].varDims.size() );
[ # # ][ # # ]
452 : :
453 : : // For a non-set variable, time should be the first dimension
454 [ # # ][ # # ]: 0 : assert( tDim == vdatas[i].varDims[0] );
[ # # ]
455 : :
456 : : // Set up readStarts and readCounts
457 [ # # ][ # # ]: 0 : vdatas[i].readStarts.resize( 3 );
458 [ # # ][ # # ]: 0 : vdatas[i].readCounts.resize( 3 );
459 : :
460 : : // First: Time
461 [ # # ][ # # ]: 0 : vdatas[i].readStarts[0] = 0; // This value is timestep dependent, will be set later
462 [ # # ][ # # ]: 0 : vdatas[i].readCounts[0] = 1;
463 : :
464 : : // Next: nVertices / nCells / nEdges
465 [ # # ]: 0 : switch( vdatas[i].entLoc )
[ # # # # ]
466 : : {
467 : : case ReadNC::ENTLOCVERT:
468 : : // Vertices
469 : : // Start from the first localGidVerts
470 : : // Actually, this will be reset later on in a loop
471 [ # # ][ # # ]: 0 : vdatas[i].readStarts[1] = localGidVerts[0] - 1;
[ # # ]
472 [ # # ][ # # ]: 0 : vdatas[i].readCounts[1] = nLocalVertices;
473 : 0 : range = &verts;
474 : 0 : break;
475 : : case ReadNC::ENTLOCFACE:
476 : : // Faces
477 : : // Start from the first localGidCells
478 : : // Actually, this will be reset later on in a loop
479 [ # # ][ # # ]: 0 : vdatas[i].readStarts[1] = localGidCells[0] - 1;
[ # # ]
480 [ # # ][ # # ]: 0 : vdatas[i].readCounts[1] = nLocalCells;
481 : 0 : range = &facesOwned;
482 : 0 : break;
483 : : case ReadNC::ENTLOCEDGE:
484 : : // Edges
485 : : // Start from the first localGidEdges
486 : : // Actually, this will be reset later on in a loop
487 [ # # ][ # # ]: 0 : vdatas[i].readStarts[1] = localGidEdges[0] - 1;
[ # # ]
488 [ # # ][ # # ]: 0 : vdatas[i].readCounts[1] = nLocalEdges;
489 : 0 : range = &edges;
490 : 0 : break;
491 : : default:
492 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << vdatas[i].varName );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
493 : : }
494 : :
495 : : // Finally: layers or other optional levels, it is possible that there is no
496 : : // level dimension (numLev is 0) for this non-set variable
497 [ # # ][ # # ]: 0 : if( vdatas[i].numLev < 1 ) vdatas[i].numLev = 1;
[ # # ]
498 [ # # ][ # # ]: 0 : vdatas[i].readStarts[2] = 0;
499 [ # # ][ # # ]: 0 : vdatas[i].readCounts[2] = vdatas[i].numLev;
[ # # ]
500 : :
501 : : // Get variable size
502 [ # # ]: 0 : vdatas[i].sz = 1;
503 [ # # ]: 0 : for( std::size_t idx = 0; idx != 3; idx++ )
504 [ # # ][ # # ]: 0 : vdatas[i].sz *= vdatas[i].readCounts[idx];
[ # # ]
505 : :
506 [ # # ]: 0 : for( unsigned int t = 0; t < tstep_nums.size(); t++ )
507 : : {
508 [ # # ][ # # ]: 0 : dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
[ # # ]
509 : :
510 [ # # ][ # # ]: 0 : if( tstep_nums[t] >= dimLens[tDim] )
[ # # ]
511 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] ); }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
512 : :
513 : : // Get the tag to read into
514 [ # # ][ # # ]: 0 : if( !vdatas[i].varTags[t] )
[ # # ]
515 : : {
516 [ # # ][ # # ]: 0 : rval = get_tag_to_nonset( vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev );MB_CHK_SET_ERR( rval, "Trouble getting tag for variable " << vdatas[i].varName );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
517 : : }
518 : :
519 : : // Get ptr to tag space
520 [ # # ][ # # ]: 0 : assert( 1 == range->psize() );
521 : : void* data;
522 : : int count;
523 [ # # ][ # # ]: 0 : rval = mbImpl->tag_iterate( vdatas[i].varTags[t], range->begin(), range->end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate tag for variable " << vdatas[i].varName );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
524 [ # # ][ # # ]: 0 : assert( (unsigned)count == range->size() );
525 [ # # ][ # # ]: 0 : vdatas[i].varDatas[t] = data;
526 : : }
527 : : }
528 : :
529 : 0 : return rval;
530 : : }
531 : :
532 : : #ifdef MOAB_HAVE_PNETCDF
533 : : ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset_async( std::vector< ReadNC::VarData >& vdatas,
534 : : std::vector< int >& tstep_nums )
535 : : {
536 : : bool& noEdges = _readNC->noEdges;
537 : : DebugOutput& dbgOut = _readNC->dbgOut;
538 : :
539 : : ErrorCode rval = read_ucd_variables_to_nonset_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read non-set variables" );
540 : :
541 : : // Finally, read into that space
542 : : int success;
543 : : Range* pLocalGid = NULL;
544 : :
545 : : for( unsigned int i = 0; i < vdatas.size(); i++ )
546 : : {
547 : : // Skip edge variables, if specified by the read options
548 : : if( noEdges && vdatas[i].entLoc == ReadNC::ENTLOCEDGE ) continue;
549 : :
550 : : switch( vdatas[i].entLoc )
551 : : {
552 : : case ReadNC::ENTLOCVERT:
553 : : pLocalGid = &localGidVerts;
554 : : break;
555 : : case ReadNC::ENTLOCFACE:
556 : : pLocalGid = &localGidCells;
557 : : break;
558 : : case ReadNC::ENTLOCEDGE:
559 : : pLocalGid = &localGidEdges;
560 : : break;
561 : : default:
562 : : MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << vdatas[i].varName );
563 : : }
564 : :
565 : : std::size_t sz = vdatas[i].sz;
566 : :
567 : : for( unsigned int t = 0; t < tstep_nums.size(); t++ )
568 : : {
569 : : // We will synchronize all these reads with the other processors,
570 : : // so the wait will be inside this double loop; is it too much?
571 : : size_t nb_reads = pLocalGid->psize();
572 : : std::vector< int > requests( nb_reads ), statuss( nb_reads );
573 : : size_t idxReq = 0;
574 : :
575 : : // Set readStart for each timestep along time dimension
576 : : vdatas[i].readStarts[0] = tstep_nums[t];
577 : :
578 : : switch( vdatas[i].varDataType )
579 : : {
580 : : case NC_FLOAT:
581 : : case NC_DOUBLE: {
582 : : // Read float as double
583 : : std::vector< double > tmpdoubledata( sz );
584 : :
585 : : // In the case of ucd mesh, and on multiple proc,
586 : : // we need to read as many times as subranges we have in the
587 : : // localGid range;
588 : : // basically, we have to give a different point
589 : : // for data to start, for every subrange :(
590 : : size_t indexInDoubleArray = 0;
591 : : size_t ic = 0;
592 : : for( Range::pair_iterator pair_iter = pLocalGid->pair_begin(); pair_iter != pLocalGid->pair_end();
593 : : ++pair_iter, ic++ )
594 : : {
595 : : EntityHandle starth = pair_iter->first;
596 : : EntityHandle endh = pair_iter->second; // inclusive
597 : : vdatas[i].readStarts[1] = ( NCDF_SIZE )( starth - 1 );
598 : : vdatas[i].readCounts[1] = ( NCDF_SIZE )( endh - starth + 1 );
599 : :
600 : : // Do a partial read, in each subrange
601 : : // wait outside this loop
602 : : success =
603 : : NCFUNCREQG( _vara_double )( _fileId, vdatas[i].varId, &( vdatas[i].readStarts[0] ),
604 : : &( vdatas[i].readCounts[0] ),
605 : : &( tmpdoubledata[indexInDoubleArray] ), &requests[idxReq++] );
606 : : if( success )
607 : : MB_SET_ERR( MB_FAILURE,
608 : : "Failed to read double data in a loop for variable " << vdatas[i].varName );
609 : : // We need to increment the index in double array for the
610 : : // next subrange
611 : : indexInDoubleArray += ( endh - starth + 1 ) * 1 * vdatas[i].numLev;
612 : : }
613 : : assert( ic == pLocalGid->psize() );
614 : :
615 : : success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
616 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
617 : :
618 : : void* data = vdatas[i].varDatas[t];
619 : : for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
620 : : ( (double*)data )[idx] = tmpdoubledata[idx];
621 : :
622 : : break;
623 : : }
624 : : default:
625 : : MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
626 : : }
627 : : }
628 : : }
629 : :
630 : : // Debug output, if requested
631 : : if( 1 == dbgOut.get_verbosity() )
632 : : {
633 : : dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
634 : : for( unsigned int i = 1; i < vdatas.size(); i++ )
635 : : dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
636 : : dbgOut.tprintf( 1, "\n" );
637 : : }
638 : :
639 : : return rval;
640 : : }
641 : : #else
642 : 0 : ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset( std::vector< ReadNC::VarData >& vdatas,
643 : : std::vector< int >& tstep_nums )
644 : : {
645 : 0 : bool& noEdges = _readNC->noEdges;
646 : 0 : DebugOutput& dbgOut = _readNC->dbgOut;
647 : :
648 [ # # ][ # # ]: 0 : ErrorCode rval = read_ucd_variables_to_nonset_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read non-set variables" );
[ # # ][ # # ]
[ # # ][ # # ]
649 : :
650 : : // Finally, read into that space
651 : : int success;
652 : 0 : Range* pLocalGid = NULL;
653 : :
654 [ # # ]: 0 : for( unsigned int i = 0; i < vdatas.size(); i++ )
655 : : {
656 : : // Skip edge variables, if specified by the read options
657 [ # # ][ # # ]: 0 : if( noEdges && vdatas[i].entLoc == ReadNC::ENTLOCEDGE ) continue;
[ # # ]
658 : :
659 [ # # # # ]: 0 : switch( vdatas[i].entLoc )
660 : : {
661 : : case ReadNC::ENTLOCVERT:
662 : 0 : pLocalGid = &localGidVerts;
663 : 0 : break;
664 : : case ReadNC::ENTLOCFACE:
665 : 0 : pLocalGid = &localGidCells;
666 : 0 : break;
667 : : case ReadNC::ENTLOCEDGE:
668 : 0 : pLocalGid = &localGidEdges;
669 : 0 : break;
670 : : default:
671 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << vdatas[i].varName );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
672 : : }
673 : :
674 : 0 : std::size_t sz = vdatas[i].sz;
675 : :
676 [ # # ]: 0 : for( unsigned int t = 0; t < tstep_nums.size(); t++ )
677 : : {
678 : : // Set readStart for each timestep along time dimension
679 : 0 : vdatas[i].readStarts[0] = tstep_nums[t];
680 : :
681 [ # # ]: 0 : switch( vdatas[i].varDataType )
682 : : {
683 : : case NC_FLOAT:
684 : : case NC_DOUBLE: {
685 : : // Read float as double
686 [ # # ]: 0 : std::vector< double > tmpdoubledata( sz );
687 : :
688 : : // In the case of ucd mesh, and on multiple proc,
689 : : // we need to read as many times as subranges we have in the
690 : : // localGid range;
691 : : // basically, we have to give a different point
692 : : // for data to start, for every subrange :(
693 : 0 : size_t indexInDoubleArray = 0;
694 : 0 : size_t ic = 0;
695 [ # # ][ # # ]: 0 : for( Range::pair_iterator pair_iter = pLocalGid->pair_begin(); pair_iter != pLocalGid->pair_end();
[ # # ][ # # ]
[ # # ]
696 : : ++pair_iter, ic++ )
697 : : {
698 [ # # ]: 0 : EntityHandle starth = pair_iter->first;
699 [ # # ]: 0 : EntityHandle endh = pair_iter->second; // Inclusive
700 [ # # ][ # # ]: 0 : vdatas[i].readStarts[1] = ( NCDF_SIZE )( starth - 1 );
701 [ # # ][ # # ]: 0 : vdatas[i].readCounts[1] = ( NCDF_SIZE )( endh - starth + 1 );
702 : :
703 [ # # ][ # # ]: 0 : success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &( vdatas[i].readStarts[0] ),
[ # # ]
704 [ # # ][ # # ]: 0 : &( vdatas[i].readCounts[0] ),
705 [ # # ][ # # ]: 0 : &( tmpdoubledata[indexInDoubleArray] ) );
706 [ # # ]: 0 : if( success )
707 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE,
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
708 : : "Failed to read double data in a loop for variable " << vdatas[i].varName );
709 : : // We need to increment the index in double array for the
710 : : // next subrange
711 [ # # ]: 0 : indexInDoubleArray += ( endh - starth + 1 ) * 1 * vdatas[i].numLev;
712 : : }
713 [ # # ][ # # ]: 0 : assert( ic == pLocalGid->psize() );
714 : :
715 [ # # ][ # # ]: 0 : void* data = vdatas[i].varDatas[t];
716 [ # # ]: 0 : for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
717 [ # # ]: 0 : ( (double*)data )[idx] = tmpdoubledata[idx];
718 : :
719 [ # # ]: 0 : break;
720 : : }
721 : : default:
722 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
723 : : }
724 : : }
725 : : }
726 : :
727 : : // Debug output, if requested
728 [ # # ]: 0 : if( 1 == dbgOut.get_verbosity() )
729 : : {
730 [ # # ][ # # ]: 0 : dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
731 [ # # ]: 0 : for( unsigned int i = 1; i < vdatas.size(); i++ )
732 : 0 : dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
733 : 0 : dbgOut.tprintf( 1, "\n" );
734 : : }
735 : :
736 : 0 : return rval;
737 : : }
738 : : #endif
739 : :
740 : : #ifdef MOAB_HAVE_MPI
741 : 0 : ErrorCode NCHelperGCRM::redistribute_local_cells( int start_cell_idx )
742 : : {
743 : : // If possible, apply Zoltan partition
744 : : #ifdef MOAB_HAVE_ZOLTAN
745 : : if( _readNC->partMethod == ScdParData::RCBZOLTAN )
746 : : {
747 : : // Read lat/lon coordinates of cell centers, then convert spherical to
748 : : // Cartesian, and use them as input to Zoltan partition
749 : : int xCellVarId;
750 : : int success = NCFUNC( inq_varid )( _fileId, "grid_center_lat", &xCellVarId );
751 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of grid_center_lat" );
752 : : std::vector< double > xCell( nLocalCells );
753 : : NCDF_SIZE read_start = static_cast< NCDF_SIZE >( start_cell_idx - 1 );
754 : : NCDF_SIZE read_count = static_cast< NCDF_SIZE >( nLocalCells );
755 : : success = NCFUNCAG( _vara_double )( _fileId, xCellVarId, &read_start, &read_count, &xCell[0] );
756 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read grid_center_lat data" );
757 : :
758 : : // Read y coordinates of cell centers
759 : : int yCellVarId;
760 : : success = NCFUNC( inq_varid )( _fileId, "grid_center_lon", &yCellVarId );
761 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of grid_center_lon" );
762 : : std::vector< double > yCell( nLocalCells );
763 : : success = NCFUNCAG( _vara_double )( _fileId, yCellVarId, &read_start, &read_count, &yCell[0] );
764 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read grid_center_lon data" );
765 : :
766 : : // Convert lon/lat/rad to x/y/z
767 : : std::vector< double > zCell( nLocalCells );
768 : : double rad = 8000.0; // This is just a approximate radius
769 : : for( int i = 0; i < nLocalCells; i++ )
770 : : {
771 : : double cosphi = cos( xCell[i] );
772 : : double zmult = sin( xCell[i] );
773 : : double xmult = cosphi * cos( yCell[i] );
774 : : double ymult = cosphi * sin( yCell[i] );
775 : : xCell[i] = rad * xmult;
776 : : yCell[i] = rad * ymult;
777 : : zCell[i] = rad * zmult;
778 : : }
779 : :
780 : : // Zoltan partition using RCB; maybe more studies would be good, as to which partition
781 : : // is better
782 : : Interface*& mbImpl = _readNC->mbImpl;
783 : : DebugOutput& dbgOut = _readNC->dbgOut;
784 : : ZoltanPartitioner* mbZTool = new ZoltanPartitioner( mbImpl, false, 0, NULL );
785 : : ErrorCode rval = mbZTool->repartition( xCell, yCell, zCell, start_cell_idx, "RCB", localGidCells );MB_CHK_SET_ERR( rval, "Error in Zoltan partitioning" );
786 : : delete mbZTool;
787 : :
788 : : dbgOut.tprintf( 1, "After Zoltan partitioning, localGidCells.psize() = %d\n", (int)localGidCells.psize() );
789 : : dbgOut.tprintf( 1, " localGidCells.size() = %d\n", (int)localGidCells.size() );
790 : :
791 : : // This is important: local cells are now redistributed, so nLocalCells might be different!
792 : : nLocalCells = localGidCells.size();
793 : :
794 : : return MB_SUCCESS;
795 : : }
796 : : #endif
797 : :
798 : : // By default, apply trivial partition
799 : 0 : localGidCells.insert( start_cell_idx, start_cell_idx + nLocalCells - 1 );
800 : :
801 : 0 : return MB_SUCCESS;
802 : : }
803 : : #endif
804 : :
805 : 0 : ErrorCode NCHelperGCRM::create_local_vertices( const std::vector< int >& vertices_on_local_cells,
806 : : EntityHandle& start_vertex )
807 : : {
808 : 0 : Interface*& mbImpl = _readNC->mbImpl;
809 : 0 : Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
810 : 0 : const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
811 : 0 : DebugOutput& dbgOut = _readNC->dbgOut;
812 : 0 : std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
813 : :
814 : : // Make a copy of vertices_on_local_cells for sorting (keep original one to set cell
815 : : // connectivity later)
816 [ # # ]: 0 : std::vector< int > vertices_on_local_cells_sorted( vertices_on_local_cells );
817 [ # # ]: 0 : std::sort( vertices_on_local_cells_sorted.begin(), vertices_on_local_cells_sorted.end() );
818 : : std::copy( vertices_on_local_cells_sorted.rbegin(), vertices_on_local_cells_sorted.rend(),
819 [ # # ][ # # ]: 0 : range_inserter( localGidVerts ) );
820 [ # # ]: 0 : nLocalVertices = localGidVerts.size();
821 : :
822 [ # # ][ # # ]: 0 : dbgOut.tprintf( 1, " localGidVerts.psize() = %d\n", (int)localGidVerts.psize() );
823 [ # # ][ # # ]: 0 : dbgOut.tprintf( 1, " localGidVerts.size() = %d\n", (int)localGidVerts.size() );
824 : :
825 : : // Create local vertices
826 [ # # ]: 0 : std::vector< double* > arrays;
827 : : ErrorCode rval =
828 : : _readNC->readMeshIface->get_node_coords( 3, nLocalVertices, 0, start_vertex, arrays,
829 : : // Might have to create gather mesh later
830 [ # # ][ # # ]: 0 : ( createGatherSet ? nLocalVertices + nVertices : nLocalVertices ) );MB_CHK_SET_ERR( rval, "Failed to create local vertices" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
831 : :
832 : : // Add local vertices to current file set
833 [ # # ]: 0 : Range local_verts_range( start_vertex, start_vertex + nLocalVertices - 1 );
834 [ # # ][ # # ]: 0 : rval = _readNC->mbImpl->add_entities( _fileSet, local_verts_range );MB_CHK_SET_ERR( rval, "Failed to add local vertices to current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
835 : :
836 : : // Get ptr to GID memory for local vertices
837 : 0 : int count = 0;
838 : 0 : void* data = NULL;
839 [ # # ][ # # ]: 0 : 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" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
840 [ # # ]: 0 : assert( count == nLocalVertices );
841 : 0 : int* gid_data = (int*)data;
842 [ # # ][ # # ]: 0 : std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
[ # # ]
843 : :
844 : : // Duplicate GID data, which will be used to resolve sharing
845 [ # # ]: 0 : if( mpFileIdTag )
846 : : {
847 [ # # ][ # # ]: 0 : 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" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
848 [ # # ]: 0 : assert( count == nLocalVertices );
849 : 0 : int bytes_per_tag = 4;
850 [ # # ][ # # ]: 0 : rval = mbImpl->tag_get_bytes( *mpFileIdTag, bytes_per_tag );MB_CHK_SET_ERR( rval, "can't get number of bytes for file id tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
851 [ # # ]: 0 : if( 4 == bytes_per_tag )
852 : : {
853 : 0 : gid_data = (int*)data;
854 [ # # ][ # # ]: 0 : std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
[ # # ]
855 : : }
856 [ # # ]: 0 : else if( 8 == bytes_per_tag )
857 : : { // Should be a handle tag on 64 bit machine?
858 : 0 : long* handle_tag_data = (long*)data;
859 [ # # ][ # # ]: 0 : std::copy( localGidVerts.begin(), localGidVerts.end(), handle_tag_data );
[ # # ]
860 : : }
861 : : }
862 : :
863 : : #ifdef MOAB_HAVE_PNETCDF
864 : : size_t nb_reads = localGidVerts.psize();
865 : : std::vector< int > requests( nb_reads );
866 : : std::vector< int > statuss( nb_reads );
867 : : size_t idxReq = 0;
868 : : #endif
869 : :
870 : : // Store lev values in levVals
871 [ # # ]: 0 : std::map< std::string, ReadNC::VarData >::iterator vmit;
872 [ # # ][ # # ]: 0 : if( ( vmit = varInfo.find( "layers" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # #
# # # # #
# ]
873 : : {
874 [ # # ][ # # ]: 0 : rval = read_coordinate( "layers", 0, nLevels - 1, levVals );MB_CHK_SET_ERR( rval, "Trouble reading 'layers' variable" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
875 : : }
876 [ # # ][ # # ]: 0 : else if( ( vmit = varInfo.find( "interfaces" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # #
# # # # #
# ]
877 : : {
878 [ # # ][ # # ]: 0 : rval = read_coordinate( "interfaces", 0, nLevels - 1, levVals );MB_CHK_SET_ERR( rval, "Trouble reading 'interfaces' variable" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
879 : : }
880 : : else
881 : : {
882 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Couldn't find 'layers' or 'interfaces' variable" );
[ # # ][ # # ]
[ # # ]
883 : : }
884 : :
885 : : // Decide whether down is positive
886 : 0 : char posval[10] = { 0 };
887 [ # # ][ # # ]: 0 : int success = NCFUNC( get_att_text )( _fileId, ( *vmit ).second.varId, "positive", posval );
888 [ # # ][ # # ]: 0 : if( 0 == success && !strncmp( posval, "down", 4 ) )
889 : : {
890 [ # # ][ # # ]: 0 : for( std::vector< double >::iterator dvit = levVals.begin(); dvit != levVals.end(); ++dvit )
[ # # ]
891 [ # # ]: 0 : ( *dvit ) *= -1.0;
892 : : }
893 : :
894 : : // Read x coordinates for local vertices
895 [ # # ]: 0 : double* xptr = arrays[0];
896 : : int xVertexVarId;
897 [ # # ]: 0 : success = NCFUNC( inq_varid )( _fileId, "grid_corner_lon", &xVertexVarId );
898 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of grid_corner_lon" );
[ # # ][ # # ]
[ # # ][ # # ]
899 : 0 : size_t indexInArray = 0;
900 [ # # ][ # # ]: 0 : for( Range::pair_iterator pair_iter = localGidVerts.pair_begin(); pair_iter != localGidVerts.pair_end();
[ # # ][ # # ]
[ # # ]
901 : : ++pair_iter )
902 : : {
903 [ # # ]: 0 : EntityHandle starth = pair_iter->first;
904 [ # # ]: 0 : EntityHandle endh = pair_iter->second;
905 : 0 : NCDF_SIZE read_start = ( NCDF_SIZE )( starth - 1 );
906 : 0 : NCDF_SIZE read_count = ( NCDF_SIZE )( endh - starth + 1 );
907 : :
908 : : // Do a partial read in each subrange
909 : : #ifdef MOAB_HAVE_PNETCDF
910 : : success = NCFUNCREQG( _vara_double )( _fileId, xVertexVarId, &read_start, &read_count, &( xptr[indexInArray] ),
911 : : &requests[idxReq++] );
912 : : #else
913 [ # # ]: 0 : success = NCFUNCAG( _vara_double )( _fileId, xVertexVarId, &read_start, &read_count, &( xptr[indexInArray] ) );
914 : : #endif
915 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read grid_corner_lon data in a loop" );
[ # # ][ # # ]
[ # # ][ # # ]
916 : :
917 : : // Increment the index for next subrange
918 : 0 : indexInArray += ( endh - starth + 1 );
919 : : }
920 : :
921 : : #ifdef MOAB_HAVE_PNETCDF
922 : : // Wait outside the loop
923 : : success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
924 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
925 : : #endif
926 : :
927 : : // Read y coordinates for local vertices
928 [ # # ]: 0 : double* yptr = arrays[1];
929 : : int yVertexVarId;
930 [ # # ]: 0 : success = NCFUNC( inq_varid )( _fileId, "grid_corner_lat", &yVertexVarId );
931 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of grid_corner_lat" );
[ # # ][ # # ]
[ # # ][ # # ]
932 : : #ifdef MOAB_HAVE_PNETCDF
933 : : idxReq = 0;
934 : : #endif
935 : 0 : indexInArray = 0;
936 [ # # ][ # # ]: 0 : for( Range::pair_iterator pair_iter = localGidVerts.pair_begin(); pair_iter != localGidVerts.pair_end();
[ # # ][ # # ]
[ # # ]
937 : : ++pair_iter )
938 : : {
939 [ # # ]: 0 : EntityHandle starth = pair_iter->first;
940 [ # # ]: 0 : EntityHandle endh = pair_iter->second;
941 : 0 : NCDF_SIZE read_start = ( NCDF_SIZE )( starth - 1 );
942 : 0 : NCDF_SIZE read_count = ( NCDF_SIZE )( endh - starth + 1 );
943 : :
944 : : // Do a partial read in each subrange
945 : : #ifdef MOAB_HAVE_PNETCDF
946 : : success = NCFUNCREQG( _vara_double )( _fileId, yVertexVarId, &read_start, &read_count, &( yptr[indexInArray] ),
947 : : &requests[idxReq++] );
948 : : #else
949 [ # # ]: 0 : success = NCFUNCAG( _vara_double )( _fileId, yVertexVarId, &read_start, &read_count, &( yptr[indexInArray] ) );
950 : : #endif
951 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read grid_corner_lat data in a loop" );
[ # # ][ # # ]
[ # # ][ # # ]
952 : :
953 : : // Increment the index for next subrange
954 : 0 : indexInArray += ( endh - starth + 1 );
955 : : }
956 : :
957 : : #ifdef MOAB_HAVE_PNETCDF
958 : : // Wait outside the loop
959 : : success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
960 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
961 : : #endif
962 : :
963 : : // Convert lon/lat/rad to x/y/z
964 [ # # ]: 0 : double* zptr = arrays[2];
965 [ # # ]: 0 : double rad = 8000.0 + levVals[0];
966 [ # # ]: 0 : for( int i = 0; i < nLocalVertices; i++ )
967 : : {
968 : 0 : double cosphi = cos( yptr[i] );
969 : 0 : double zmult = sin( yptr[i] );
970 : 0 : double xmult = cosphi * cos( xptr[i] );
971 : 0 : double ymult = cosphi * sin( xptr[i] );
972 : 0 : xptr[i] = rad * xmult;
973 : 0 : yptr[i] = rad * ymult;
974 : 0 : zptr[i] = rad * zmult;
975 : : }
976 : :
977 : 0 : return MB_SUCCESS;
978 : : }
979 : :
980 : 0 : ErrorCode NCHelperGCRM::create_local_edges( EntityHandle start_vertex )
981 : : {
982 : 0 : Interface*& mbImpl = _readNC->mbImpl;
983 : 0 : Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
984 : 0 : DebugOutput& dbgOut = _readNC->dbgOut;
985 : :
986 : : // Read edges on each local cell, to get localGidEdges
987 : : int edgesOnCellVarId;
988 [ # # ]: 0 : int success = NCFUNC( inq_varid )( _fileId, "cell_edges", &edgesOnCellVarId );
989 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of cell_edges" );
[ # # ][ # # ]
[ # # ][ # # ]
990 : :
991 [ # # ]: 0 : std::vector< int > edges_on_local_cells( nLocalCells * EDGES_PER_CELL );
992 [ # # ]: 0 : dbgOut.tprintf( 1, " edges_on_local_cells.size() = %d\n", (int)edges_on_local_cells.size() );
993 : :
994 : : #ifdef MOAB_HAVE_PNETCDF
995 : : size_t nb_reads = localGidCells.psize();
996 : : std::vector< int > requests( nb_reads );
997 : : std::vector< int > statuss( nb_reads );
998 : : size_t idxReq = 0;
999 : : #endif
1000 : 0 : size_t indexInArray = 0;
1001 [ # # ][ # # ]: 0 : for( Range::pair_iterator pair_iter = localGidCells.pair_begin(); pair_iter != localGidCells.pair_end();
[ # # ][ # # ]
[ # # ]
1002 : : ++pair_iter )
1003 : : {
1004 [ # # ]: 0 : EntityHandle starth = pair_iter->first;
1005 [ # # ]: 0 : EntityHandle endh = pair_iter->second;
1006 [ # # ]: 0 : dbgOut.tprintf( 1, " starth = %d\n", (int)starth );
1007 [ # # ]: 0 : dbgOut.tprintf( 1, " endh = %d\n", (int)endh );
1008 : 0 : NCDF_SIZE read_starts[2] = { static_cast< NCDF_SIZE >( starth - 1 ), 0 };
1009 : 0 : NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( endh - starth + 1 ),
1010 : 0 : static_cast< NCDF_SIZE >( EDGES_PER_CELL ) };
1011 : :
1012 : : // Do a partial read in each subrange
1013 : : #ifdef MOAB_HAVE_PNETCDF
1014 : : success = NCFUNCREQG( _vara_int )( _fileId, edgesOnCellVarId, read_starts, read_counts,
1015 : : &( edges_on_local_cells[indexInArray] ), &requests[idxReq++] );
1016 : : #else
1017 : : success = NCFUNCAG( _vara_int )( _fileId, edgesOnCellVarId, read_starts, read_counts,
1018 [ # # ][ # # ]: 0 : &( edges_on_local_cells[indexInArray] ) );
1019 : : #endif
1020 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read cell_edges data in a loop" );
[ # # ][ # # ]
[ # # ][ # # ]
1021 : :
1022 : : // Increment the index for next subrange
1023 : 0 : indexInArray += ( endh - starth + 1 ) * EDGES_PER_CELL;
1024 : : }
1025 : :
1026 : : #ifdef MOAB_HAVE_PNETCDF
1027 : : // Wait outside the loop
1028 : : success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
1029 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
1030 : : #endif
1031 : :
1032 : : // GCRM is 0 based, convert edge indices from 0 to 1 based
1033 [ # # ]: 0 : for( std::size_t idx = 0; idx < edges_on_local_cells.size(); idx++ )
1034 [ # # ]: 0 : edges_on_local_cells[idx] += 1;
1035 : :
1036 : : // Collect local edges
1037 [ # # ]: 0 : std::sort( edges_on_local_cells.begin(), edges_on_local_cells.end() );
1038 [ # # ][ # # ]: 0 : std::copy( edges_on_local_cells.rbegin(), edges_on_local_cells.rend(), range_inserter( localGidEdges ) );
1039 [ # # ]: 0 : nLocalEdges = localGidEdges.size();
1040 : :
1041 [ # # ][ # # ]: 0 : dbgOut.tprintf( 1, " localGidEdges.psize() = %d\n", (int)localGidEdges.psize() );
1042 [ # # ][ # # ]: 0 : dbgOut.tprintf( 1, " localGidEdges.size() = %d\n", (int)localGidEdges.size() );
1043 : :
1044 : : // Create local edges
1045 : : EntityHandle start_edge;
1046 : 0 : EntityHandle* conn_arr_edges = NULL;
1047 : : ErrorCode rval =
1048 : : _readNC->readMeshIface->get_element_connect( nLocalEdges, 2, MBEDGE, 0, start_edge, conn_arr_edges,
1049 : : // Might have to create gather mesh later
1050 [ # # ][ # # ]: 0 : ( createGatherSet ? nLocalEdges + nEdges : nLocalEdges ) );MB_CHK_SET_ERR( rval, "Failed to create local edges" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1051 : :
1052 : : // Add local edges to current file set
1053 [ # # ]: 0 : Range local_edges_range( start_edge, start_edge + nLocalEdges - 1 );
1054 [ # # ][ # # ]: 0 : rval = _readNC->mbImpl->add_entities( _fileSet, local_edges_range );MB_CHK_SET_ERR( rval, "Failed to add local edges to current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1055 : :
1056 : : // Get ptr to GID memory for edges
1057 : 0 : int count = 0;
1058 : 0 : void* data = NULL;
1059 [ # # ][ # # ]: 0 : rval = mbImpl->tag_iterate( mGlobalIdTag, local_edges_range.begin(), local_edges_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate global id tag on local edges" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1060 [ # # ]: 0 : assert( count == nLocalEdges );
1061 : 0 : int* gid_data = (int*)data;
1062 [ # # ][ # # ]: 0 : std::copy( localGidEdges.begin(), localGidEdges.end(), gid_data );
[ # # ]
1063 : :
1064 : : int verticesOnEdgeVarId;
1065 : :
1066 : : // Read vertices on each local edge, to get edge connectivity
1067 [ # # ]: 0 : success = NCFUNC( inq_varid )( _fileId, "edge_corners", &verticesOnEdgeVarId );
1068 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of edge_corners" );
[ # # ][ # # ]
[ # # ][ # # ]
1069 : : // Utilize the memory storage pointed by conn_arr_edges
1070 : 0 : int* vertices_on_local_edges = (int*)conn_arr_edges;
1071 : : #ifdef MOAB_HAVE_PNETCDF
1072 : : nb_reads = localGidEdges.psize();
1073 : : requests.resize( nb_reads );
1074 : : statuss.resize( nb_reads );
1075 : : idxReq = 0;
1076 : : #endif
1077 : 0 : indexInArray = 0;
1078 [ # # ][ # # ]: 0 : for( Range::pair_iterator pair_iter = localGidEdges.pair_begin(); pair_iter != localGidEdges.pair_end();
[ # # ][ # # ]
[ # # ]
1079 : : ++pair_iter )
1080 : : {
1081 [ # # ]: 0 : EntityHandle starth = pair_iter->first;
1082 [ # # ]: 0 : EntityHandle endh = pair_iter->second;
1083 : 0 : NCDF_SIZE read_starts[2] = { static_cast< NCDF_SIZE >( starth - 1 ), 0 };
1084 : 0 : NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( endh - starth + 1 ), 2 };
1085 : :
1086 : : // Do a partial read in each subrange
1087 : : #ifdef MOAB_HAVE_PNETCDF
1088 : : success = NCFUNCREQG( _vara_int )( _fileId, verticesOnEdgeVarId, read_starts, read_counts,
1089 : : &( vertices_on_local_edges[indexInArray] ), &requests[idxReq++] );
1090 : : #else
1091 : : success = NCFUNCAG( _vara_int )( _fileId, verticesOnEdgeVarId, read_starts, read_counts,
1092 [ # # ]: 0 : &( vertices_on_local_edges[indexInArray] ) );
1093 : : #endif
1094 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read edge_corners data in a loop" );
[ # # ][ # # ]
[ # # ][ # # ]
1095 : :
1096 : : // Increment the index for next subrange
1097 : 0 : indexInArray += ( endh - starth + 1 ) * 2;
1098 : : }
1099 : :
1100 : : #ifdef MOAB_HAVE_PNETCDF
1101 : : // Wait outside the loop
1102 : : success = NCFUNC( wait_all )( _fileId, requests.size(), &requests[0], &statuss[0] );
1103 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
1104 : : #endif
1105 : :
1106 : : // Populate connectivity data for local edges
1107 : : // Convert in-place from int (stored in the first half) to EntityHandle
1108 : : // Reading backward is the trick
1109 [ # # ]: 0 : for( int edge_vert = nLocalEdges * 2 - 1; edge_vert >= 0; edge_vert-- )
1110 : : {
1111 : : // Note, indices stored in vertices_on_local_edges are 0 based
1112 : 0 : int global_vert_idx = vertices_on_local_edges[edge_vert] + 1; // Global vertex index, 1 based
1113 [ # # ]: 0 : int local_vert_idx = localGidVerts.index( global_vert_idx ); // Local vertex index, 0 based
1114 [ # # ]: 0 : assert( local_vert_idx != -1 );
1115 : 0 : conn_arr_edges[edge_vert] = start_vertex + local_vert_idx;
1116 : : }
1117 : :
1118 : 0 : return MB_SUCCESS;
1119 : : }
1120 : :
1121 : 0 : ErrorCode NCHelperGCRM::create_padded_local_cells( const std::vector< int >& vertices_on_local_cells,
1122 : : EntityHandle start_vertex, Range& faces )
1123 : : {
1124 : 0 : Interface*& mbImpl = _readNC->mbImpl;
1125 : 0 : Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
1126 : :
1127 : : // Create cells
1128 : : EntityHandle start_element;
1129 : 0 : EntityHandle* conn_arr_local_cells = NULL;
1130 : : ErrorCode rval =
1131 : : _readNC->readMeshIface->get_element_connect( nLocalCells, EDGES_PER_CELL, MBPOLYGON, 0, start_element,
1132 : : conn_arr_local_cells,
1133 : : // Might have to create gather mesh later
1134 [ # # ][ # # ]: 0 : ( createGatherSet ? nLocalCells + nCells : nLocalCells ) );MB_CHK_SET_ERR( rval, "Failed to create local cells" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1135 [ # # ]: 0 : faces.insert( start_element, start_element + nLocalCells - 1 );
1136 : :
1137 : : // Add local cells to current file set
1138 [ # # ]: 0 : Range local_cells_range( start_element, start_element + nLocalCells - 1 );
1139 [ # # ][ # # ]: 0 : rval = _readNC->mbImpl->add_entities( _fileSet, local_cells_range );MB_CHK_SET_ERR( rval, "Failed to add local cells to current file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1140 : :
1141 : : // Get ptr to GID memory for local cells
1142 : 0 : int count = 0;
1143 : 0 : void* data = NULL;
1144 [ # # ][ # # ]: 0 : 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" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1145 [ # # ]: 0 : assert( count == nLocalCells );
1146 : 0 : int* gid_data = (int*)data;
1147 [ # # ][ # # ]: 0 : std::copy( localGidCells.begin(), localGidCells.end(), gid_data );
[ # # ]
1148 : :
1149 : : // Set connectivity array with proper local vertices handles
1150 : : // vertices_on_local_cells array was already corrected to have
1151 : : // the last vertices repeated for pentagons, e.g. 122345 => 123455
1152 [ # # ]: 0 : for( int local_cell_idx = 0; local_cell_idx < nLocalCells; local_cell_idx++ )
1153 : : {
1154 [ # # ]: 0 : for( int i = 0; i < EDGES_PER_CELL; i++ )
1155 : : {
1156 : : // Note, indices stored in vertices_on_local_cells are 1 based
1157 : : EntityHandle global_vert_idx =
1158 [ # # ]: 0 : vertices_on_local_cells[local_cell_idx * EDGES_PER_CELL + i]; // Global vertex index, 1 based
1159 [ # # ]: 0 : int local_vert_idx = localGidVerts.index( global_vert_idx ); // Local vertex index, 0 based
1160 [ # # ]: 0 : assert( local_vert_idx != -1 );
1161 : 0 : conn_arr_local_cells[local_cell_idx * EDGES_PER_CELL + i] = start_vertex + local_vert_idx;
1162 : : }
1163 : : }
1164 : :
1165 : 0 : return MB_SUCCESS;
1166 : : }
1167 : :
1168 : 0 : ErrorCode NCHelperGCRM::create_gather_set_vertices( EntityHandle gather_set, EntityHandle& gather_set_start_vertex )
1169 : : {
1170 : 0 : Interface*& mbImpl = _readNC->mbImpl;
1171 : 0 : Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
1172 : 0 : const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
1173 : :
1174 : : // Create gather set vertices
1175 [ # # ]: 0 : std::vector< double* > arrays;
1176 : : // Don't need to specify allocation number here, because we know enough vertices were created
1177 : : // before
1178 [ # # ][ # # ]: 0 : ErrorCode rval = _readNC->readMeshIface->get_node_coords( 3, nVertices, 0, gather_set_start_vertex, arrays );MB_CHK_SET_ERR( rval, "Failed to create gather set vertices" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1179 : :
1180 : : // Add vertices to the gather set
1181 [ # # ]: 0 : Range gather_set_verts_range( gather_set_start_vertex, gather_set_start_vertex + nVertices - 1 );
1182 [ # # ][ # # ]: 0 : rval = mbImpl->add_entities( gather_set, gather_set_verts_range );MB_CHK_SET_ERR( rval, "Failed to add vertices to the gather set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1183 : :
1184 : : // Read x coordinates for gather set vertices
1185 [ # # ]: 0 : double* xptr = arrays[0];
1186 : : int xVertexVarId;
1187 [ # # ]: 0 : int success = NCFUNC( inq_varid )( _fileId, "grid_corner_lon", &xVertexVarId );
1188 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of grid_corner_lon" );
[ # # ][ # # ]
[ # # ][ # # ]
1189 : 0 : NCDF_SIZE read_start = 0;
1190 : 0 : NCDF_SIZE read_count = static_cast< NCDF_SIZE >( nVertices );
1191 : : #ifdef MOAB_HAVE_PNETCDF
1192 : : // Enter independent I/O mode, since this read is only for the gather processor
1193 : : success = NCFUNC( begin_indep_data )( _fileId );
1194 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to begin independent I/O mode" );
1195 : : success = NCFUNCG( _vara_double )( _fileId, xVertexVarId, &read_start, &read_count, xptr );
1196 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read grid_corner_lon data" );
1197 : : success = NCFUNC( end_indep_data )( _fileId );
1198 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to end independent I/O mode" );
1199 : : #else
1200 [ # # ]: 0 : success = NCFUNCG( _vara_double )( _fileId, xVertexVarId, &read_start, &read_count, xptr );
1201 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read grid_corner_lon data" );
[ # # ][ # # ]
[ # # ][ # # ]
1202 : : #endif
1203 : :
1204 : : // Read y coordinates for gather set vertices
1205 [ # # ]: 0 : double* yptr = arrays[1];
1206 : : int yVertexVarId;
1207 [ # # ]: 0 : success = NCFUNC( inq_varid )( _fileId, "grid_corner_lat", &yVertexVarId );
1208 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of grid_corner_lat" );
[ # # ][ # # ]
[ # # ][ # # ]
1209 : : #ifdef MOAB_HAVE_PNETCDF
1210 : : // Enter independent I/O mode, since this read is only for the gather processor
1211 : : success = NCFUNC( begin_indep_data )( _fileId );
1212 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to begin independent I/O mode" );
1213 : : success = NCFUNCG( _vara_double )( _fileId, yVertexVarId, &read_start, &read_count, yptr );
1214 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read grid_corner_lat data" );
1215 : : success = NCFUNC( end_indep_data )( _fileId );
1216 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to end independent I/O mode" );
1217 : : #else
1218 [ # # ]: 0 : success = NCFUNCG( _vara_double )( _fileId, yVertexVarId, &read_start, &read_count, yptr );
1219 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read grid_corner_lat data" );
[ # # ][ # # ]
[ # # ][ # # ]
1220 : : #endif
1221 : :
1222 : : // Convert lon/lat/rad to x/y/z
1223 [ # # ]: 0 : double* zptr = arrays[2];
1224 [ # # ]: 0 : double rad = 8000.0 + levVals[0];
1225 [ # # ]: 0 : for( int i = 0; i < nVertices; i++ )
1226 : : {
1227 : 0 : double cosphi = cos( yptr[i] );
1228 : 0 : double zmult = sin( yptr[i] );
1229 : 0 : double xmult = cosphi * cos( xptr[i] );
1230 : 0 : double ymult = cosphi * sin( xptr[i] );
1231 : 0 : xptr[i] = rad * xmult;
1232 : 0 : yptr[i] = rad * ymult;
1233 : 0 : zptr[i] = rad * zmult;
1234 : : }
1235 : :
1236 : : // Get ptr to GID memory for gather set vertices
1237 : 0 : int count = 0;
1238 : 0 : void* data = NULL;
1239 : : rval =
1240 [ # # ][ # # ]: 0 : mbImpl->tag_iterate( mGlobalIdTag, gather_set_verts_range.begin(), gather_set_verts_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate global id tag on gather set vertices" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1241 [ # # ]: 0 : assert( count == nVertices );
1242 : 0 : int* gid_data = (int*)data;
1243 [ # # ]: 0 : for( int j = 1; j <= nVertices; j++ )
1244 : 0 : gid_data[j - 1] = j;
1245 : :
1246 : : // Set the file id tag too, it should be bigger something not interfering with global id
1247 [ # # ]: 0 : if( mpFileIdTag )
1248 : : {
1249 : : rval = mbImpl->tag_iterate( *mpFileIdTag, gather_set_verts_range.begin(), gather_set_verts_range.end(), count,
1250 [ # # ][ # # ]: 0 : data );MB_CHK_SET_ERR( rval, "Failed to iterate file id tag on gather set vertices" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1251 [ # # ]: 0 : assert( count == nVertices );
1252 : 0 : int bytes_per_tag = 4;
1253 [ # # ][ # # ]: 0 : rval = mbImpl->tag_get_bytes( *mpFileIdTag, bytes_per_tag );MB_CHK_SET_ERR( rval, "Can't get number of bytes for file id tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1254 [ # # ]: 0 : if( 4 == bytes_per_tag )
1255 : : {
1256 : 0 : gid_data = (int*)data;
1257 [ # # ]: 0 : for( int j = 1; j <= nVertices; j++ )
1258 : 0 : gid_data[j - 1] = nVertices + j; // Bigger than global id tag
1259 : : }
1260 [ # # ]: 0 : else if( 8 == bytes_per_tag )
1261 : : { // Should be a handle tag on 64 bit machine?
1262 : 0 : long* handle_tag_data = (long*)data;
1263 [ # # ]: 0 : for( int j = 1; j <= nVertices; j++ )
1264 : 0 : handle_tag_data[j - 1] = nVertices + j; // Bigger than global id tag
1265 : : }
1266 : : }
1267 : :
1268 : 0 : return MB_SUCCESS;
1269 : : }
1270 : :
1271 : 0 : ErrorCode NCHelperGCRM::create_gather_set_edges( EntityHandle gather_set, EntityHandle gather_set_start_vertex )
1272 : : {
1273 : 0 : Interface*& mbImpl = _readNC->mbImpl;
1274 : :
1275 : : // Create gather set edges
1276 : : EntityHandle start_edge;
1277 : 0 : EntityHandle* conn_arr_gather_set_edges = NULL;
1278 : : // Don't need to specify allocation number here, because we know enough edges were created
1279 : : // before
1280 : : ErrorCode rval =
1281 [ # # ][ # # ]: 0 : _readNC->readMeshIface->get_element_connect( nEdges, 2, MBEDGE, 0, start_edge, conn_arr_gather_set_edges );MB_CHK_SET_ERR( rval, "Failed to create gather set edges" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1282 : :
1283 : : // Add edges to the gather set
1284 [ # # ]: 0 : Range gather_set_edges_range( start_edge, start_edge + nEdges - 1 );
1285 [ # # ][ # # ]: 0 : rval = mbImpl->add_entities( gather_set, gather_set_edges_range );MB_CHK_SET_ERR( rval, "Failed to add edges to the gather set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1286 : :
1287 : : // Read vertices on each edge
1288 : : int verticesOnEdgeVarId;
1289 [ # # ]: 0 : int success = NCFUNC( inq_varid )( _fileId, "edge_corners", &verticesOnEdgeVarId );
1290 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of edge_corners" );
[ # # ][ # # ]
[ # # ][ # # ]
1291 : : // Utilize the memory storage pointed by conn_arr_gather_set_edges
1292 : 0 : int* vertices_on_gather_set_edges = (int*)conn_arr_gather_set_edges;
1293 : 0 : NCDF_SIZE read_starts[2] = { 0, 0 };
1294 : 0 : NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( nEdges ), 2 };
1295 : : #ifdef MOAB_HAVE_PNETCDF
1296 : : // Enter independent I/O mode, since this read is only for the gather processor
1297 : : success = NCFUNC( begin_indep_data )( _fileId );
1298 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to begin independent I/O mode" );
1299 : : success =
1300 : : NCFUNCG( _vara_int )( _fileId, verticesOnEdgeVarId, read_starts, read_counts, vertices_on_gather_set_edges );
1301 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read edge_corners data" );
1302 : : success = NCFUNC( end_indep_data )( _fileId );
1303 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to end independent I/O mode" );
1304 : : #else
1305 : : success =
1306 [ # # ]: 0 : NCFUNCG( _vara_int )( _fileId, verticesOnEdgeVarId, read_starts, read_counts, vertices_on_gather_set_edges );
1307 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read edge_corners data" );
[ # # ][ # # ]
[ # # ][ # # ]
1308 : : #endif
1309 : :
1310 : : // Populate connectivity data for gather set edges
1311 : : // Convert in-place from int (stored in the first half) to EntityHandle
1312 : : // Reading backward is the trick
1313 [ # # ]: 0 : for( int edge_vert = nEdges * 2 - 1; edge_vert >= 0; edge_vert-- )
1314 : : {
1315 : : // Note, indices stored in vertices_on_gather_set_edges are 0 based
1316 : 0 : int gather_set_vert_idx = vertices_on_gather_set_edges[edge_vert]; // Global vertex index, 0 based
1317 : : // Connectivity array is shifted by where the gather set vertices start
1318 : 0 : conn_arr_gather_set_edges[edge_vert] = gather_set_start_vertex + gather_set_vert_idx;
1319 : : }
1320 : :
1321 : 0 : return MB_SUCCESS;
1322 : : }
1323 : :
1324 : 0 : ErrorCode NCHelperGCRM::create_padded_gather_set_cells( EntityHandle gather_set, EntityHandle gather_set_start_vertex )
1325 : : {
1326 : 0 : Interface*& mbImpl = _readNC->mbImpl;
1327 : :
1328 : : // Create gather set cells
1329 : : EntityHandle start_element;
1330 : 0 : EntityHandle* conn_arr_gather_set_cells = NULL;
1331 : : // Don't need to specify allocation number here, because we know enough cells were created
1332 : : // before
1333 : : ErrorCode rval = _readNC->readMeshIface->get_element_connect( nCells, EDGES_PER_CELL, MBPOLYGON, 0, start_element,
1334 [ # # ][ # # ]: 0 : conn_arr_gather_set_cells );MB_CHK_SET_ERR( rval, "Failed to create gather set cells" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1335 : :
1336 : : // Add cells to the gather set
1337 [ # # ]: 0 : Range gather_set_cells_range( start_element, start_element + nCells - 1 );
1338 [ # # ][ # # ]: 0 : rval = mbImpl->add_entities( gather_set, gather_set_cells_range );MB_CHK_SET_ERR( rval, "Failed to add cells to the gather set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1339 : :
1340 : : // Read vertices on each gather set cell (connectivity)
1341 : : int verticesOnCellVarId;
1342 [ # # ]: 0 : int success = NCFUNC( inq_varid )( _fileId, "cell_corners", &verticesOnCellVarId );
1343 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of cell_corners" );
[ # # ][ # # ]
[ # # ][ # # ]
1344 : : // Utilize the memory storage pointed by conn_arr_gather_set_cells
1345 : 0 : int* vertices_on_gather_set_cells = (int*)conn_arr_gather_set_cells;
1346 : 0 : NCDF_SIZE read_starts[2] = { 0, 0 };
1347 : 0 : NCDF_SIZE read_counts[2] = { static_cast< NCDF_SIZE >( nCells ), static_cast< NCDF_SIZE >( EDGES_PER_CELL ) };
1348 : : #ifdef MOAB_HAVE_PNETCDF
1349 : : // Enter independent I/O mode, since this read is only for the gather processor
1350 : : success = NCFUNC( begin_indep_data )( _fileId );
1351 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to begin independent I/O mode" );
1352 : : success =
1353 : : NCFUNCG( _vara_int )( _fileId, verticesOnCellVarId, read_starts, read_counts, vertices_on_gather_set_cells );
1354 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read cell_corners data" );
1355 : : success = NCFUNC( end_indep_data )( _fileId );
1356 : : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to end independent I/O mode" );
1357 : : #else
1358 : : success =
1359 [ # # ]: 0 : NCFUNCG( _vara_int )( _fileId, verticesOnCellVarId, read_starts, read_counts, vertices_on_gather_set_cells );
1360 [ # # ][ # # ]: 0 : if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read cell_corners data" );
[ # # ][ # # ]
[ # # ][ # # ]
1361 : : #endif
1362 : :
1363 : : // Correct gather set cell vertices array in the same way as local cell vertices array
1364 : : // Pentagons as hexagons should have a connectivity like 123455 and not 122345
1365 [ # # ]: 0 : for( int gather_set_cell_idx = 0; gather_set_cell_idx < nCells; gather_set_cell_idx++ )
1366 : : {
1367 : 0 : int* pvertex = vertices_on_gather_set_cells + gather_set_cell_idx * EDGES_PER_CELL;
1368 [ # # ]: 0 : for( int k = 0; k < EDGES_PER_CELL - 2; k++ )
1369 : : {
1370 [ # # ]: 0 : if( *( pvertex + k ) == *( pvertex + k + 1 ) )
1371 : : {
1372 : : // Shift the connectivity
1373 [ # # ]: 0 : for( int kk = k + 1; kk < EDGES_PER_CELL - 1; kk++ )
1374 : 0 : *( pvertex + kk ) = *( pvertex + kk + 1 );
1375 : : // No need to try next k
1376 : 0 : break;
1377 : : }
1378 : : }
1379 : : }
1380 : :
1381 : : // Populate connectivity data for gather set cells
1382 : : // Convert in-place from int (stored in the first half) to EntityHandle
1383 : : // Reading backward is the trick
1384 [ # # ]: 0 : for( int cell_vert = nCells * EDGES_PER_CELL - 1; cell_vert >= 0; cell_vert-- )
1385 : : {
1386 : : // Note, indices stored in vertices_on_gather_set_cells are 0 based
1387 : 0 : int gather_set_vert_idx = vertices_on_gather_set_cells[cell_vert]; // Global vertex index, 0 based
1388 : : // Connectivity array is shifted by where the gather set vertices start
1389 : 0 : conn_arr_gather_set_cells[cell_vert] = gather_set_start_vertex + gather_set_vert_idx;
1390 : : }
1391 : :
1392 : 0 : return MB_SUCCESS;
1393 : : }
1394 : :
1395 [ + - ][ + - ]: 228 : } // namespace moab
|
