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