 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "NCHelper.hpp"
2 #include "NCHelperEuler.hpp"
3 #include "NCHelperFV.hpp"
4 #include "NCHelperDomain.hpp"
5 #include "NCHelperScrip.hpp"
6 #include "NCHelperHOMME.hpp"
7 #include "NCHelperMPAS.hpp"
8 #include "NCHelperGCRM.hpp"
9 #include "NCHelperESMF.hpp"
10
11 #include <sstream>
12
13 #include "MBTagConventions.hpp"
14
15 #ifdef WIN32
16 #ifdef size_t
17 #undef size_t
18 #endif
19 #endif
20
21 namespace moab
22 {
23
24 ReadNC::NCFormatType NCHelper::get_nc_format( ReadNC* readNC, int fileId )
25 {
26 // Check if CF convention is being followed
27 bool is_CF = false;
28
29 std::map< std::string, ReadNC::AttData >& globalAtts = readNC->globalAtts;
30 std::map< std::string, ReadNC::AttData >::iterator attIt = globalAtts.find( "conventions" );
31 if( attIt == globalAtts.end() ) attIt = globalAtts.find( "Conventions" );
32
33 if( attIt != globalAtts.end() )
34 {
35 unsigned int sz = attIt->second.attLen;
36 std::string att_data;
37 att_data.resize( sz + 1 );
38 att_data[sz] = '\000';
39 int success =
40 NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
41 if( 0 == success && att_data.find( "CF" ) != std::string::npos ) is_CF = true;
42 }
43
44 // Depending on the format, update FileOptions as well
45 if( NCHelperMPAS::can_read_file( readNC ) )
46 return ReadNC::NC_FORMAT_MPAS;
47 else if( NCHelperScrip::can_read_file( readNC, fileId ) )
48 return ReadNC::NC_FORMAT_SCRIP;
49 else if( NCHelperESMF::can_read_file( readNC ) )
50 return ReadNC::NC_FORMAT_ESMF;
51 else if( NCHelperDomain::can_read_file( readNC, fileId ) ) // && is_CF )
52 return ReadNC::NC_FORMAT_DOMAIN;
53 else if( NCHelperHOMME::can_read_file( readNC, fileId ) )
54 return ReadNC::NC_FORMAT_HOMME;
55 else if( NCHelperEuler::can_read_file( readNC, fileId ) && is_CF )
56 return ReadNC::NC_FORMAT_EULER;
57 else if( NCHelperGCRM::can_read_file( readNC ) )
58 return ReadNC::NC_FORMAT_GCRM;
59 else if( NCHelperFV::can_read_file( readNC, fileId ) && is_CF )
60 return ReadNC::NC_FORMAT_FV;
61 else // Unknown NetCDF grid (will fill this in later for POP, CICE and CLM)
62 return ReadNC::NC_FORMAT_UNKNOWN_TYPE;
63 }
64
65 //! Get appropriate file read options depending on the format of the NC file
66 std::string NCHelper::get_default_ncformat_options( ReadNC::NCFormatType format )
67 {
68 switch( format )
69 {
70 case moab::ReadNC::NC_FORMAT_GCRM: // GCRM format reader
71 case moab::ReadNC::NC_FORMAT_MPAS: // MPAS format reader
72 return "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;"
73 "PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
74 case moab::ReadNC::NC_FORMAT_SCRIP: // SCRIP format reader
75 return "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;";
76 case moab::ReadNC::NC_FORMAT_ESMF: // ESMF unstructured format reader
77 return "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;PARALLEL_RESOLVE_SHARED_ENTS;VARIABLE=;";
78 case moab::ReadNC::NC_FORMAT_DOMAIN: // Climate Domain reader
79 return "PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;";
80 case moab::ReadNC::NC_FORMAT_HOMME: // HOMME format reader
81 return "PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL;PARALLEL_RESOLVE_SHARED_ENTS;";
82 case moab::ReadNC::NC_FORMAT_EULER: // Euler format reader
83 case moab::ReadNC::NC_FORMAT_FV: // FV climate format reader
84 return "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
85 "PARTITION_METHOD=SQIJ;VARIABLE=;";
86 default:
87 return "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;PARALLEL_RESOLVE_SHARED_ENTS;";
88 }
89 }
90
91 NCHelper* NCHelper::get_nc_helper( ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet )
92 {
93 ReadNC::NCFormatType nctype = NCHelper::get_nc_format( readNC, fileId );
94
95 switch( nctype )
96 {
97 case ReadNC::NC_FORMAT_MPAS: // MPAS format reader
98 return new( std::nothrow ) NCHelperMPAS( readNC, fileId, opts, fileSet );
99 case ReadNC::NC_FORMAT_SCRIP: // SCRIP format reader
100 // SCRIP can be CF or non CF, if it comes from MPAS :)
101 return new( std::nothrow ) NCHelperScrip( readNC, fileId, opts, fileSet );
102 case ReadNC::NC_FORMAT_ESMF: // ESMF unstructured format reader
103 return new( std::nothrow ) NCHelperESMF( readNC, fileId, opts, fileSet );
104 case ReadNC::NC_FORMAT_DOMAIN: // Climate Domain reader
105 return new( std::nothrow ) NCHelperDomain( readNC, fileId, opts, fileSet );
106 case ReadNC::NC_FORMAT_HOMME: // HOMME format reader
107 // For a HOMME connectivity file, there might be no CF convention
108 return new( std::nothrow ) NCHelperHOMME( readNC, fileId, opts, fileSet );
109 case ReadNC::NC_FORMAT_GCRM: // GCRM format reader
110 return new( std::nothrow ) NCHelperGCRM( readNC, fileId, opts, fileSet );
111 case ReadNC::NC_FORMAT_EULER: // Euler format reader
112 return new( std::nothrow ) NCHelperEuler( readNC, fileId, opts, fileSet );
113 case ReadNC::NC_FORMAT_FV: // FV climate format reader
114 return new( std::nothrow ) NCHelperFV( readNC, fileId, opts, fileSet );
115 default: // Unknown NetCDF grid (will fill this in later for POP, CICE and CLM)
116 return nullptr;
117 }
118 }
119
120 ErrorCode NCHelper::create_conventional_tags( const std::vector< int >& tstep_nums )
121 {
122 Interface*& mbImpl = _readNC->mbImpl;
123 std::vector< std::string >& dimNames = _readNC->dimNames;
124 std::vector< int >& dimLens = _readNC->dimLens;
125 std::map< std::string, ReadNC::AttData >& globalAtts = _readNC->globalAtts;
126 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
127 DebugOutput& dbgOut = _readNC->dbgOut;
128 int& partMethod = _readNC->partMethod;
129 ScdInterface* scdi = _readNC->scdi;
130
131 ErrorCode rval;
132 std::string tag_name;
133
134 // <__NUM_DIMS>
135 Tag numDimsTag = 0;
136 tag_name = "__NUM_DIMS";
137 int numDims = dimNames.size();
138 rval = mbImpl->tag_get_handle( tag_name.c_str(), 1, MB_TYPE_INTEGER, numDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
139 rval = mbImpl->tag_set_data( numDimsTag, &_fileSet, 1, &numDims );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
140 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
141
142 // <__NUM_VARS>
143 Tag numVarsTag = 0;
144 tag_name = "__NUM_VARS";
145 int numVars = varInfo.size();
146 rval = mbImpl->tag_get_handle( tag_name.c_str(), 1, MB_TYPE_INTEGER, numVarsTag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
147 rval = mbImpl->tag_set_data( numVarsTag, &_fileSet, 1, &numVars );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
148 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
149
150 // <__DIM_NAMES>
151 Tag dimNamesTag = 0;
152 tag_name = "__DIM_NAMES";
153 std::string dimnames;
154 unsigned int dimNamesSz = dimNames.size();
155 for( unsigned int i = 0; i != dimNamesSz; i++ )
156 {
157 dimnames.append( dimNames[i] );
158 dimnames.push_back( '\0' );
159 }
160 int dimnamesSz = dimnames.size();
161 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, dimNamesTag,
162 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
163 const void* ptr = dimnames.c_str();
164 rval = mbImpl->tag_set_by_ptr( dimNamesTag, &_fileSet, 1, &ptr, &dimnamesSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
165 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
166
167 // <__DIM_LENS>
168 Tag dimLensTag = 0;
169 tag_name = "__DIM_LENS";
170 int dimLensSz = dimLens.size();
171 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, dimLensTag,
172 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
173 ptr = &( dimLens[0] );
174 rval = mbImpl->tag_set_by_ptr( dimLensTag, &_fileSet, 1, &ptr, &dimLensSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
175 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
176
177 // <__VAR_NAMES>
178 Tag varNamesTag = 0;
179 tag_name = "__VAR_NAMES";
180 std::string varnames;
181 std::map< std::string, ReadNC::VarData >::iterator mapIter;
182 for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
183 {
184 varnames.append( mapIter->first );
185 varnames.push_back( '\0' );
186 }
187 int varnamesSz = varnames.size();
188 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varNamesTag,
189 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
190 ptr = varnames.c_str();
191 rval = mbImpl->tag_set_by_ptr( varNamesTag, &_fileSet, 1, &ptr, &varnamesSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
192 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
193
194 // __<dim_name>_LOC_MINMAX (for time)
195 for( unsigned int i = 0; i != dimNamesSz; i++ )
196 {
197 if( dimNames[i] == "time" || dimNames[i] == "Time" || dimNames[i] == "t" )
198 {
199 // some files might have Time dimension as 0, it will not appear in any
200 // variables, so skip it
201 if( nTimeSteps == 0 ) continue;
202 std::stringstream ss_tag_name;
203 ss_tag_name << "__" << dimNames[i] << "_LOC_MINMAX";
204 tag_name = ss_tag_name.str();
205 Tag tagh = 0;
206 std::vector< int > val( 2, 0 );
207 val[0] = 0;
208 val[1] = nTimeSteps - 1;
209 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 );
210 rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
211 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
212 }
213 }
214
215 // __<dim_name>_LOC_VALS (for time)
216 for( unsigned int i = 0; i != dimNamesSz; i++ )
217 {
218 if( dimNames[i] == "time" || dimNames[i] == "Time" || dimNames[i] == "t" )
219 {
220 std::vector< int > val;
221 if( !tstep_nums.empty() )
222 val = tstep_nums;
223 else
224 {
225 // some files might have Time dimension as 0, it will not appear in any
226 // variables, so skip it
227 if( tVals.empty() ) continue;
228 val.resize( tVals.size() );
229 for( unsigned int j = 0; j != tVals.size(); j++ )
230 val[j] = j;
231 }
232 Tag tagh = 0;
233 std::stringstream ss_tag_name;
234 ss_tag_name << "__" << dimNames[i] << "_LOC_VALS";
235 tag_name = ss_tag_name.str();
236 rval = mbImpl->tag_get_handle( tag_name.c_str(), val.size(), MB_TYPE_INTEGER, tagh,
237 MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
238 rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
239 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
240 }
241 }
242
243 // __<var_name>_DIMS
244 for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
245 {
246 Tag varNamesDimsTag = 0;
247 std::stringstream ss_tag_name;
248 ss_tag_name << "__" << mapIter->first << "_DIMS";
249 tag_name = ss_tag_name.str();
250 unsigned int varDimSz = varInfo[mapIter->first].varDims.size();
251 if( varDimSz == 0 ) continue;
252 std::vector< Tag > varDimTags( varDimSz );
253 for( unsigned int i = 0; i != varDimSz; i++ )
254 {
255 Tag tmptag = 0;
256 std::string tmptagname = dimNames[varInfo[mapIter->first].varDims[i]];
257 rval = mbImpl->tag_get_handle( tmptagname.c_str(), 0, MB_TYPE_OPAQUE, tmptag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting tag " << tmptagname );
258 varDimTags[i] = tmptag;
259 }
260 // rval = mbImpl->tag_get_handle(tag_name.c_str(), varDimSz, MB_TYPE_HANDLE,
261 // varNamesDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT); We should not use MB_TYPE_HANDLE for Tag
262 // here. Tag is a pointer, which is 4 bytes on 32 bit machines and 8 bytes on 64 bit
263 // machines. Normally, entity handle is 8 bytes on 64 bit machines, but it can also be
264 // configured to 4 bytes.
265 rval = mbImpl->tag_get_handle( tag_name.c_str(), varDimSz * sizeof( Tag ), MB_TYPE_OPAQUE, varNamesDimsTag,
266 MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
267 rval = mbImpl->tag_set_data( varNamesDimsTag, &_fileSet, 1, &( varDimTags[0] ) );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
268 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
269 }
270
271 // <PARTITION_METHOD>
272 Tag part_tag = scdi->part_method_tag();
273 if( !part_tag ) MB_SET_ERR( MB_FAILURE, "Trouble getting PARTITION_METHOD tag" );
274 rval = mbImpl->tag_set_data( part_tag, &_fileSet, 1, &partMethod );MB_CHK_SET_ERR( rval, "Trouble setting data to PARTITION_METHOD tag" );
275 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
276
277 // <__GLOBAL_ATTRIBS>
278 tag_name = "__GLOBAL_ATTRIBS";
279 Tag globalAttTag = 0;
280 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, globalAttTag,
281 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
282 std::string gattVal;
283 std::vector< int > gattLen;
284 rval = create_attrib_string( globalAtts, gattVal, gattLen );MB_CHK_SET_ERR( rval, "Trouble creating global attribute string" );
285 const void* gattptr = gattVal.c_str();
286 int globalAttSz = gattVal.size();
287 rval = mbImpl->tag_set_by_ptr( globalAttTag, &_fileSet, 1, &gattptr, &globalAttSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
288 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
289
290 // <__GLOBAL_ATTRIBS_LEN>
291 tag_name = "__GLOBAL_ATTRIBS_LEN";
292 Tag globalAttLenTag = 0;
293 if( gattLen.size() == 0 ) gattLen.push_back( 0 );
294 rval = mbImpl->tag_get_handle( tag_name.c_str(), gattLen.size(), MB_TYPE_INTEGER, globalAttLenTag,
295 MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
296 rval = mbImpl->tag_set_data( globalAttLenTag, &_fileSet, 1, &gattLen[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
297 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
298
299 // __<var_name>_ATTRIBS and __<var_name>_ATTRIBS_LEN
300 for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
301 {
302 std::stringstream ssTagName;
303 ssTagName << "__" << mapIter->first << "_ATTRIBS";
304 tag_name = ssTagName.str();
305 Tag varAttTag = 0;
306 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varAttTag,
307 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
308
309 std::string varAttVal;
310 std::vector< int > varAttLen;
311 if( mapIter->second.numAtts < 1 )
312 {
313 if( dummyVarNames.find( mapIter->first ) != dummyVarNames.end() )
314 {
315 // This variable is a dummy coordinate variable
316 varAttVal = "DUMMY_VAR";
317 }
318 else
319 {
320 // This variable has no attributes
321 varAttVal = "NO_ATTRIBS";
322 }
323 }
324 else
325 {
326 rval = create_attrib_string( mapIter->second.varAtts, varAttVal, varAttLen );MB_CHK_SET_ERR( rval, "Trouble creating attribute string for variable " << mapIter->first );
327 }
328 const void* varAttPtr = varAttVal.c_str();
329 int varAttSz = varAttVal.size();
330 if( 0 == varAttSz ) varAttSz = 1;
331 rval = mbImpl->tag_set_by_ptr( varAttTag, &_fileSet, 1, &varAttPtr, &varAttSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
332 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
333
334 ssTagName << "_LEN";
335 tag_name = ssTagName.str();
336 Tag varAttLenTag = 0;
337 if( 0 == varAttLen.size() ) varAttLen.push_back( 0 );
338 rval = mbImpl->tag_get_handle( tag_name.c_str(), varAttLen.size(), MB_TYPE_INTEGER, varAttLenTag,
339 MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
340 rval = mbImpl->tag_set_data( varAttLenTag, &_fileSet, 1, &varAttLen[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
341 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
342 }
343
344 // <__VAR_NAMES_LOCATIONS>
345 tag_name = "__VAR_NAMES_LOCATIONS";
346 Tag varNamesLocsTag = 0;
347 std::vector< int > varNamesLocs( varInfo.size() );
348 rval = mbImpl->tag_get_handle( tag_name.c_str(), varNamesLocs.size(), MB_TYPE_INTEGER, varNamesLocsTag,
349 MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
350 for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
351 {
352 varNamesLocs[std::distance( varInfo.begin(), mapIter )] = mapIter->second.entLoc;
353 }
354 rval = mbImpl->tag_set_data( varNamesLocsTag, &_fileSet, 1, &varNamesLocs[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
355 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
356
357 // <__MESH_TYPE>
358 Tag meshTypeTag = 0;
359 tag_name = "__MESH_TYPE";
360 std::string meshTypeName = get_mesh_type_name();
361
362 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, meshTypeTag,
363 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
364 ptr = meshTypeName.c_str();
365 int leng = meshTypeName.size();
366 rval = mbImpl->tag_set_by_ptr( meshTypeTag, &_fileSet, 1, &ptr, &leng );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
367 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
368
369 return MB_SUCCESS;
370 }
371
372 ErrorCode NCHelper::update_time_tag_vals()
373 {
374 Interface*& mbImpl = _readNC->mbImpl;
375 std::vector< std::string >& dimNames = _readNC->dimNames;
376
377 ErrorCode rval;
378
379 // The time tag might be a dummy one (e.g. 'Time' for MPAS)
380 std::string time_tag_name = dimNames[tDim];
381 if( dummyVarNames.find( time_tag_name ) != dummyVarNames.end() ) return MB_SUCCESS;
382
383 Tag time_tag = 0;
384 const void* data = NULL;
385 int time_tag_size = 0;
386 rval = mbImpl->tag_get_handle( time_tag_name.c_str(), 0, MB_TYPE_DOUBLE, time_tag, MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble getting tag " << time_tag_name );
387 rval = mbImpl->tag_get_by_ptr( time_tag, &_fileSet, 1, &data, &time_tag_size );MB_CHK_SET_ERR( rval, "Trouble getting data of tag " << time_tag_name );
388 const double* time_tag_vals = static_cast< const double* >( data );
389
390 // Merge tVals (read from current file) to existing time tag
391 // Assume that time_tag_vals and tVals are both sorted
392 std::vector< double > merged_time_vals;
393 merged_time_vals.reserve( time_tag_size + nTimeSteps );
394 int i = 0;
395 int j = 0;
396
397 // Merge time values from time_tag_vals and tVals
398 while( i < time_tag_size && j < nTimeSteps )
399 {
400 if( time_tag_vals[i] < tVals[j] )
401 merged_time_vals.push_back( time_tag_vals[i++] );
402 else
403 merged_time_vals.push_back( tVals[j++] );
404 }
405
406 // Append remaining time values of time_tag_vals (if any)
407 while( i < time_tag_size )
408 merged_time_vals.push_back( time_tag_vals[i++] );
409
410 // Append remaining time values of tVals (if any)
411 while( j < nTimeSteps )
412 merged_time_vals.push_back( tVals[j++] );
413
414 data = &merged_time_vals[0];
415 time_tag_size = merged_time_vals.size();
416 rval = mbImpl->tag_set_by_ptr( time_tag, &_fileSet, 1, &data, &time_tag_size );MB_CHK_SET_ERR( rval, "Trouble setting data to tag " << time_tag_name );
417
418 return MB_SUCCESS;
419 }
420
421 ErrorCode NCHelper::read_variables_setup( std::vector< std::string >& var_names,
422 std::vector< int >& tstep_nums,
423 std::vector< ReadNC::VarData >& vdatas,
424 std::vector< ReadNC::VarData >& vsetdatas )
425 {
426 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
427 std::vector< std::string >& dimNames = _readNC->dimNames;
428
429 std::map< std::string, ReadNC::VarData >::iterator mit;
430
431 // If empty read them all (except ignored variables)
432 if( var_names.empty() )
433 {
434 for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
435 {
436 ReadNC::VarData vd = ( *mit ).second;
437
438 // If read all variables at once, skip ignored ones
439 if( ignoredVarNames.find( vd.varName ) != ignoredVarNames.end() ) continue;
440
441 // Coordinate variables (include dummy ones) were read to the file set by default
442 if( std::find( dimNames.begin(), dimNames.end(), vd.varName ) != dimNames.end() ) continue;
443
444 if( vd.entLoc == ReadNC::ENTLOCSET )
445 vsetdatas.push_back( vd );
446 else
447 vdatas.push_back( vd );
448 }
449 }
450 else
451 {
452 // Read specified variables (might include ignored ones)
453 for( unsigned int i = 0; i < var_names.size(); i++ )
454 {
455 mit = varInfo.find( var_names[i] );
456 if( mit != varInfo.end() )
457 {
458 ReadNC::VarData vd = ( *mit ).second;
459
460 // Upon creation of dummy coordinate variables, tag values have already been set
461 if( dummyVarNames.find( vd.varName ) != dummyVarNames.end() ) continue;
462
463 if( vd.entLoc == ReadNC::ENTLOCSET )
464 vsetdatas.push_back( vd );
465 else
466 vdatas.push_back( vd );
467 }
468 else
469 {
470 MB_SET_ERR( MB_FAILURE, "Couldn't find specified variable " << var_names[i] );
471 }
472 }
473 }
474
475 if( tstep_nums.empty() && nTimeSteps > 0 )
476 {
477 // No timesteps input, get them all
478 for( int i = 0; i < nTimeSteps; i++ )
479 tstep_nums.push_back( i );
480 }
481
482 if( !tstep_nums.empty() )
483 {
484 for( unsigned int i = 0; i < vdatas.size(); i++ )
485 {
486 vdatas[i].varTags.resize( tstep_nums.size(), 0 );
487 vdatas[i].varDatas.resize( tstep_nums.size() );
488 // NC reader assumes that non-set variables always have timesteps
489 assert( std::find( vdatas[i].varDims.begin(), vdatas[i].varDims.end(), tDim ) != vdatas[i].varDims.end() );
490 vdatas[i].has_tsteps = true;
491 }
492
493 for( unsigned int i = 0; i < vsetdatas.size(); i++ )
494 {
495 if( ( std::find( vsetdatas[i].varDims.begin(), vsetdatas[i].varDims.end(), tDim ) !=
496 vsetdatas[i].varDims.end() ) &&
497 ( vsetdatas[i].varName != dimNames[tDim] ) )
498 {
499 // Set variables with timesteps: e.g. xtime(Time) or xtime(Time, StrLen)
500 vsetdatas[i].varTags.resize( tstep_nums.size(), 0 );
501 vsetdatas[i].varDatas.resize( tstep_nums.size() );
502 vsetdatas[i].has_tsteps = true;
503 }
504 else
505 {
506 // Set variables without timesteps: no time dimension, or time itself
507 vsetdatas[i].varTags.resize( 1, 0 );
508 vsetdatas[i].varDatas.resize( 1 );
509 vsetdatas[i].has_tsteps = false;
510 }
511 }
512 }
513
514 return MB_SUCCESS;
515 }
516
517 ErrorCode NCHelper::read_variables_to_set( std::vector< ReadNC::VarData >& vdatas, std::vector< int >& tstep_nums )
518 {
519 Interface*& mbImpl = _readNC->mbImpl;
520 DebugOutput& dbgOut = _readNC->dbgOut;
521
522 ErrorCode rval = read_variables_to_set_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read set variables" );
523
524 // Finally, read into that space
525 int success;
526 for( unsigned int i = 0; i < vdatas.size(); i++ )
527 {
528 // Note, for set variables without timesteps, loop one time and then break
529 for( unsigned int t = 0; t < tstep_nums.size(); t++ )
530 {
531 void* data = vdatas[i].varDatas[t];
532
533 // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen)
534 if( vdatas[i].has_tsteps )
535 {
536 // Set readStart for each timestep along time dimension
537 vdatas[i].readStarts[0] = tstep_nums[t];
538 }
539
540 switch( vdatas[i].varDataType )
541 {
542 case NC_BYTE:
543 case NC_CHAR:
544 success = NCFUNCAG( _vara_text )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
545 &vdatas[i].readCounts[0], (char*)data );
546 if( success )
547 MB_SET_ERR( MB_FAILURE, "Failed to read byte/char data for variable " << vdatas[i].varName );
548 break;
549 case NC_SHORT:
550 case NC_INT:
551 success = NCFUNCAG( _vara_int )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
552 &vdatas[i].readCounts[0], (int*)data );
553 if( success )
554 MB_SET_ERR( MB_FAILURE, "Failed to read short/int data for variable " << vdatas[i].varName );
555 break;
556 case NC_INT64:
557 success = NCFUNCAG( _vara_long )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
558 &vdatas[i].readCounts[0], (long*)data );
559 if( success )
560 MB_SET_ERR( MB_FAILURE, "Failed to read long data for variable " << vdatas[i].varName );
561 break;
562 case NC_FLOAT:
563 case NC_DOUBLE:
564 success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
565 &vdatas[i].readCounts[0], (double*)data );
566 if( success )
567 MB_SET_ERR( MB_FAILURE, "Failed to read float/double data for variable " << vdatas[i].varName );
568 break;
569 default:
570 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
571 }
572
573 dbgOut.tprintf( 2, "Setting data for variable %s, time step %d\n", vdatas[i].varName.c_str(),
574 tstep_nums[t] );
575 rval = mbImpl->tag_set_by_ptr( vdatas[i].varTags[t], &_fileSet, 1, &data, &vdatas[i].sz );MB_CHK_SET_ERR( rval, "Trouble setting tag data for variable " << vdatas[i].varName );
576
577 // Memory pointed by pointer data can be deleted, as tag_set_by_ptr() has already copied
578 // the tag values
579 switch( vdatas[i].varDataType )
580 {
581 case NC_BYTE:
582 case NC_CHAR:
583 delete[](char*)data;
584 break;
585 case NC_SHORT:
586 case NC_INT:
587 delete[](int*)data;
588 break;
589 case NC_INT64:
590 delete[](long*)data;
591 break;
592 case NC_FLOAT:
593 case NC_DOUBLE:
594 delete[](double*)data;
595 break;
596 default:
597 break;
598 }
599 vdatas[i].varDatas[t] = NULL;
600
601 // Loop continues only for set variables with timesteps, e.g. xtime(Time) or xtime(Time,
602 // StrLen)
603 if( !vdatas[i].has_tsteps ) break;
604 }
605 }
606
607 // Debug output, if requested
608 if( 1 == dbgOut.get_verbosity() )
609 {
610 dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
611 for( unsigned int i = 1; i < vdatas.size(); i++ )
612 dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
613 dbgOut.tprintf( 1, "\n" );
614 }
615
616 return rval;
617 }
618
619 ErrorCode NCHelper::read_coordinate( const char* var_name, int lmin, int lmax, std::vector< double >& cvals )
620 {
621 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
622 std::map< std::string, ReadNC::VarData >::iterator vmit = varInfo.find( var_name );
623 if( varInfo.end() == vmit ) MB_SET_ERR( MB_FAILURE, "Couldn't find variable " << var_name );
624
625 assert( lmin >= 0 && lmax >= lmin );
626 NCDF_SIZE tstart = lmin;
627 NCDF_SIZE tcount = lmax - lmin + 1;
628 NCDF_DIFF dum_stride = 1;
629 int success;
630
631 // Check size
632 if( (std::size_t)tcount != cvals.size() ) cvals.resize( tcount );
633
634 // Check to make sure it's a float or double
635 switch( ( *vmit ).second.varDataType )
636 {
637 case NC_FLOAT:
638 case NC_DOUBLE:
639 // Read float as double
640 success =
641 NCFUNCAG( _vars_double )( _fileId, ( *vmit ).second.varId, &tstart, &tcount, &dum_stride, &cvals[0] );
642 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read float/double data for variable " << var_name );
643 break;
644 default:
645 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << var_name );
646 }
647
648 return MB_SUCCESS;
649 }
650
651 ErrorCode NCHelper::get_tag_to_set( ReadNC::VarData& var_data, int tstep_num, Tag& tagh )
652 {
653 Interface*& mbImpl = _readNC->mbImpl;
654 DebugOutput& dbgOut = _readNC->dbgOut;
655 int& tStepBase = _readNC->tStepBase;
656
657 if( tStepBase > 0 ) tstep_num += tStepBase;
658
659 std::ostringstream tag_name;
660 if( var_data.has_tsteps )
661 tag_name << var_data.varName << tstep_num;
662 else
663 tag_name << var_data.varName;
664
665 ErrorCode rval = MB_SUCCESS;
666 tagh = 0;
667 switch( var_data.varDataType )
668 {
669 case NC_BYTE:
670 case NC_CHAR:
671 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), 0, MB_TYPE_OPAQUE, tagh,
672 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
673 break;
674 case NC_SHORT:
675 case NC_INT:
676 case NC_INT64: // this is a big stretch; we should introduce LONG tag
677 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), 0, MB_TYPE_INTEGER, tagh,
678 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
679 break;
680 case NC_FLOAT:
681 case NC_DOUBLE:
682 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), 0, MB_TYPE_DOUBLE, tagh,
683 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
684 break;
685 default:
686 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << var_data.varName );
687 }
688
689 dbgOut.tprintf( 2, "Tag %s created\n", tag_name.str().c_str() );
690
691 return rval;
692 }
693
694 ErrorCode NCHelper::get_tag_to_nonset( ReadNC::VarData& var_data, int tstep_num, Tag& tagh, int num_lev )
695 {
696 Interface*& mbImpl = _readNC->mbImpl;
697 DebugOutput& dbgOut = _readNC->dbgOut;
698 int& tStepBase = _readNC->tStepBase;
699
700 if( tStepBase > 0 ) tstep_num += tStepBase;
701
702 std::ostringstream tag_name;
703 tag_name << var_data.varName << tstep_num;
704
705 ErrorCode rval = MB_SUCCESS;
706 tagh = 0;
707 switch( var_data.varDataType )
708 {
709 case NC_BYTE:
710 case NC_CHAR:
711 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), num_lev, MB_TYPE_OPAQUE, tagh,
712 MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
713 break;
714 case NC_SHORT:
715 case NC_INT:
716 case NC_INT64:
717 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), num_lev, MB_TYPE_INTEGER, tagh,
718 MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
719 break;
720 case NC_FLOAT:
721 case NC_DOUBLE:
722 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), num_lev, MB_TYPE_DOUBLE, tagh,
723 MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
724 break;
725 default:
726 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << var_data.varName );
727 }
728
729 dbgOut.tprintf( 2, "Tag %s created\n", tag_name.str().c_str() );
730
731 return rval;
732 }
733
734 ErrorCode NCHelper::create_attrib_string( const std::map< std::string, ReadNC::AttData >& attMap,
735 std::string& attVal,
736 std::vector< int >& attLen )
737 {
738 int success;
739 std::stringstream ssAtt;
740 unsigned int sz = 0;
741 std::map< std::string, ReadNC::AttData >::const_iterator attIt = attMap.begin();
742 for( ; attIt != attMap.end(); ++attIt )
743 {
744 ssAtt << attIt->second.attName;
745 ssAtt << '\0';
746 void* attData = NULL;
747 switch( attIt->second.attDataType )
748 {
749 case NC_BYTE:
750 case NC_CHAR:
751 sz = attIt->second.attLen;
752 attData = (char*)malloc( sz );
753 success = NCFUNC( get_att_text )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
754 (char*)attData );
755 if( success )
756 MB_SET_ERR( MB_FAILURE, "Failed to read byte/char data for attribute " << attIt->second.attName );
757 ssAtt << "char;";
758 break;
759 case NC_SHORT:
760 sz = attIt->second.attLen * sizeof( short );
761 attData = (short*)malloc( sz );
762 success = NCFUNC( get_att_short )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
763 (short*)attData );
764 if( success )
765 MB_SET_ERR( MB_FAILURE, "Failed to read short data for attribute " << attIt->second.attName );
766 ssAtt << "short;";
767 break;
768 case NC_INT:
769 sz = attIt->second.attLen * sizeof( int );
770 attData = (int*)malloc( sz );
771 success = NCFUNC( get_att_int )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
772 (int*)attData );
773 if( success )
774 MB_SET_ERR( MB_FAILURE, "Failed to read int data for attribute " << attIt->second.attName );
775 ssAtt << "int;";
776 break;
777 case NC_INT64: // be careful here
778 sz = attIt->second.attLen * sizeof( long );
779 attData = (long*)malloc( sz );
780 success = NCFUNC( get_att_long )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
781 (long*)attData );
782 if( success )
783 MB_SET_ERR( MB_FAILURE, "Failed to read int data for attribute " << attIt->second.attName );
784 ssAtt << "long;";
785 break;
786 case NC_FLOAT:
787 sz = attIt->second.attLen * sizeof( float );
788 attData = (float*)malloc( sz );
789 success = NCFUNC( get_att_float )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
790 (float*)attData );
791 if( success )
792 MB_SET_ERR( MB_FAILURE, "Failed to read float data for attribute " << attIt->second.attName );
793 ssAtt << "float;";
794 break;
795 case NC_DOUBLE:
796 sz = attIt->second.attLen * sizeof( double );
797 attData = (double*)malloc( sz );
798 success = NCFUNC( get_att_double )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
799 (double*)attData );
800 if( success )
801 MB_SET_ERR( MB_FAILURE, "Failed to read double data for attribute " << attIt->second.attName );
802 ssAtt << "double;";
803 break;
804 default:
805 MB_SET_ERR( MB_FAILURE, "Unexpected data type for attribute " << attIt->second.attName );
806 }
807 char* tmpc = (char*)attData;
808 for( unsigned int counter = 0; counter != sz; ++counter )
809 ssAtt << tmpc[counter];
810 free( attData );
811 ssAtt << ';';
812 attLen.push_back( ssAtt.str().size() - 1 );
813 }
814 attVal = ssAtt.str();
815
816 return MB_SUCCESS;
817 }
818
819 ErrorCode NCHelper::create_dummy_variables()
820 {
821 Interface*& mbImpl = _readNC->mbImpl;
822 std::vector< std::string >& dimNames = _readNC->dimNames;
823 std::vector< int >& dimLens = _readNC->dimLens;
824 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
825 DebugOutput& dbgOut = _readNC->dbgOut;
826
827 // Hack: look at all dimensions, and see if we have one that does not appear in the list of
828 // varInfo names Right now, candidates are from unstructured meshes, such as ncol (HOMME) and
829 // nCells (MPAS) For each of them, create a dummy coordinate variable with a sparse tag to store
830 // the dimension length
831 for( unsigned int i = 0; i < dimNames.size(); i++ )
832 {
833 // If there is a variable with this dimension name, skip
834 if( varInfo.find( dimNames[i] ) != varInfo.end() ) continue;
835
836 // Create a dummy coordinate variable
837 int sizeTotalVar = varInfo.size();
838 std::string var_name( dimNames[i] );
839 ReadNC::VarData& data = varInfo[var_name];
840 data.varName = var_name;
841 data.varId = sizeTotalVar;
842 data.varTags.resize( 1, 0 );
843 data.varDataType = NC_INT;
844 data.varDims.resize( 1 );
845 data.varDims[0] = (int)i;
846 data.numAtts = 0;
847 data.entLoc = ReadNC::ENTLOCSET;
848 dummyVarNames.insert( var_name );
849 dbgOut.tprintf( 2, "Dummy coordinate variable created for dimension %s\n", var_name.c_str() );
850
851 // Create a corresponding sparse tag
852 Tag tagh;
853 ErrorCode rval = mbImpl->tag_get_handle( var_name.c_str(), 0, MB_TYPE_INTEGER, tagh,
854 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag for dummy coordinate variable " << var_name );
855
856 // Tag value is the dimension length
857 const void* ptr = &dimLens[i];
858 // Tag size is 1
859 int size = 1;
860 rval = mbImpl->tag_set_by_ptr( tagh, &_fileSet, 1, &ptr, &size );MB_CHK_SET_ERR( rval, "Trouble setting tag data for dummy coordinate variable " << var_name );
861
862 dbgOut.tprintf( 2, "Sparse tag created for dimension %s\n", var_name.c_str() );
863 }
864
865 return MB_SUCCESS;
866 }
867
868 ErrorCode NCHelper::read_variables_to_set_allocate( std::vector< ReadNC::VarData >& vdatas,
869 std::vector< int >& tstep_nums )
870 {
871 std::vector< int >& dimLens = _readNC->dimLens;
872 DebugOutput& dbgOut = _readNC->dbgOut;
873
874 ErrorCode rval = MB_SUCCESS;
875
876 for( unsigned int i = 0; i < vdatas.size(); i++ )
877 {
878 // Set up readStarts and readCounts
879 if( vdatas[i].has_tsteps )
880 {
881 // First: time
882 vdatas[i].readStarts.push_back( 0 ); // This value is timestep dependent, will be set later
883 vdatas[i].readCounts.push_back( 1 );
884
885 // Next: other dimensions
886 for( unsigned int idx = 1; idx != vdatas[i].varDims.size(); idx++ )
887 {
888 vdatas[i].readStarts.push_back( 0 );
889 vdatas[i].readCounts.push_back( dimLens[vdatas[i].varDims[idx]] );
890 }
891 }
892 else
893 {
894 if( vdatas[i].varDims.empty() )
895 {
896 // Scalar variable
897 vdatas[i].readStarts.push_back( 0 );
898 vdatas[i].readCounts.push_back( 1 );
899 }
900 else
901 {
902 for( unsigned int idx = 0; idx != vdatas[i].varDims.size(); idx++ )
903 {
904 vdatas[i].readStarts.push_back( 0 );
905 vdatas[i].readCounts.push_back( dimLens[vdatas[i].varDims[idx]] );
906 }
907 }
908 }
909
910 // Get variable size
911 vdatas[i].sz = 1;
912 for( std::size_t idx = 0; idx != vdatas[i].readCounts.size(); idx++ )
913 vdatas[i].sz *= vdatas[i].readCounts[idx];
914
915 // Note, for set variables without timesteps, loop one time and then break
916 for( unsigned int t = 0; t < tstep_nums.size(); t++ )
917 {
918 dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
919
920 if( tstep_nums[t] >= dimLens[tDim] )
921 {
922 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] );
923 }
924
925 // Get the tag to read into
926 if( !vdatas[i].varTags[t] )
927 {
928 rval = get_tag_to_set( vdatas[i], tstep_nums[t], vdatas[i].varTags[t] );MB_CHK_SET_ERR( rval, "Trouble getting tag to set variable " << vdatas[i].varName );
929 }
930
931 switch( vdatas[i].varDataType )
932 {
933 case NC_BYTE:
934 case NC_CHAR:
935 vdatas[i].varDatas[t] = new char[vdatas[i].sz];
936 break;
937 case NC_SHORT:
938 case NC_INT:
939 vdatas[i].varDatas[t] = new int[vdatas[i].sz];
940 break;
941 case NC_INT64:
942 vdatas[i].varDatas[t] = new long[vdatas[i].sz];
943 break;
944 case NC_FLOAT:
945 case NC_DOUBLE:
946 vdatas[i].varDatas[t] = new double[vdatas[i].sz];
947 break;
948 default:
949 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
950 }
951
952 // Loop continues only for set variables with timesteps, e.g. xtime(Time) or xtime(Time,
953 // StrLen)
954 if( !vdatas[i].has_tsteps ) break;
955 }
956 }
957
958 return rval;
959 }
960
961 ErrorCode ScdNCHelper::check_existing_mesh()
962 {
963 Interface*& mbImpl = _readNC->mbImpl;
964
965 // Get the number of vertices
966 int num_verts;
967 ErrorCode rval = mbImpl->get_number_entities_by_dimension( _fileSet, 0, num_verts );MB_CHK_SET_ERR( rval, "Trouble getting number of vertices" );
968
969 /*
970 // Check against parameters
971 // When ghosting is used, this check might fail (to be updated later)
972 if (num_verts > 0) {
973 int expected_verts = (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1) * (-1 == lDims[2] ?
974 1 : lDims[5] - lDims[2] + 1); if (num_verts != expected_verts) { MB_SET_ERR(MB_FAILURE, "Number
975 of vertices doesn't match");
976 }
977 }
978 */
979
980 // Check the number of elements too
981 int num_elems;
982 rval = mbImpl->get_number_entities_by_dimension( _fileSet, ( -1 == lCDims[2] ? 2 : 3 ), num_elems );MB_CHK_SET_ERR( rval, "Trouble getting number of elements" );
983
984 /*
985 // Check against parameters
986 // When ghosting is used, this check might fail (to be updated later)
987 if (num_elems > 0) {
988 int expected_elems = (lCDims[3] - lCDims[0] + 1) * (lCDims[4] - lCDims[1] + 1) * (-1 ==
989 lCDims[2] ? 1 : (lCDims[5] - lCDims[2] + 1)); if (num_elems != expected_elems) {
990 MB_SET_ERR(MB_FAILURE, "Number of elements doesn't match");
991 }
992 }
993 */
994
995 return MB_SUCCESS;
996 }
997
998 ErrorCode ScdNCHelper::create_mesh( Range& faces )
999 {
1000 Interface*& mbImpl = _readNC->mbImpl;
1001 Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
1002 const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
1003 DebugOutput& dbgOut = _readNC->dbgOut;
1004 ScdInterface* scdi = _readNC->scdi;
1005 ScdParData& parData = _readNC->parData;
1006
1007 Range tmp_range;
1008 ScdBox* scd_box;
1009
1010 ErrorCode rval =
1011 scdi->construct_box( HomCoord( lDims[0], lDims[1], lDims[2], 1 ), HomCoord( lDims[3], lDims[4], lDims[5], 1 ),
1012 NULL, 0, scd_box, locallyPeriodic, &parData, true );MB_CHK_SET_ERR( rval, "Trouble creating scd vertex sequence" );
1013
1014 // Add verts to tmp_range first, so we can duplicate global ids in vertex ids
1015 tmp_range.insert( scd_box->start_vertex(), scd_box->start_vertex() + scd_box->num_vertices() - 1 );
1016
1017 if( mpFileIdTag )
1018 {
1019 int count;
1020 void* data;
1021 rval = mbImpl->tag_iterate( *mpFileIdTag, tmp_range.begin(), tmp_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate file ID tag on local vertices" );
1022 assert( count == scd_box->num_vertices() );
1023 int* fid_data = (int*)data;
1024 rval = mbImpl->tag_iterate( mGlobalIdTag, tmp_range.begin(), tmp_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate global ID tag on local vertices" );
1025 assert( count == scd_box->num_vertices() );
1026 int* gid_data = (int*)data;
1027 for( int i = 0; i < count; i++ )
1028 fid_data[i] = gid_data[i];
1029 }
1030
1031 // Then add box set and elements to the range, then to the file set
1032 tmp_range.insert( scd_box->start_element(), scd_box->start_element() + scd_box->num_elements() - 1 );
1033 tmp_range.insert( scd_box->box_set() );
1034 rval = mbImpl->add_entities( _fileSet, tmp_range );MB_CHK_SET_ERR( rval, "Couldn't add new vertices to current file set" );
1035
1036 dbgOut.tprintf( 1, "scdbox %d quads, %d vertices\n", scd_box->num_elements(), scd_box->num_vertices() );
1037
1038 // Set the vertex coordinates
1039 double *xc, *yc, *zc;
1040 rval = scd_box->get_coordinate_arrays( xc, yc, zc );MB_CHK_SET_ERR( rval, "Couldn't get vertex coordinate arrays" );
1041
1042 int i, j, k, il, jl, kl;
1043 int dil = lDims[3] - lDims[0] + 1;
1044 int djl = lDims[4] - lDims[1] + 1;
1045 assert( dil == (int)ilVals.size() && djl == (int)jlVals.size() &&
1046 ( -1 == lDims[2] || lDims[5] - lDims[2] + 1 == (int)levVals.size() ) );
1047
1048 for( kl = lDims[2]; kl <= lDims[5]; kl++ )
1049 {
1050 k = kl - lDims[2];
1051 for( jl = lDims[1]; jl <= lDims[4]; jl++ )
1052 {
1053 j = jl - lDims[1];
1054 for( il = lDims[0]; il <= lDims[3]; il++ )
1055 {
1056 i = il - lDims[0];
1057 unsigned int pos = i + j * dil + k * dil * djl;
1058 xc[pos] = ilVals[i];
1059 yc[pos] = jlVals[j];
1060 zc[pos] = ( -1 == lDims[2] ? 0.0 : levVals[k] );
1061 }
1062 }
1063 }
1064
1065 #ifndef NDEBUG
1066 int num_verts =
1067 ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) * ( -1 == lDims[2] ? 1 : lDims[5] - lDims[2] + 1 );
1068 std::vector< int > gids( num_verts );
1069 Range verts( scd_box->start_vertex(), scd_box->start_vertex() + scd_box->num_vertices() - 1 );
1070 rval = mbImpl->tag_get_data( mGlobalIdTag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting local gid values of vertices" );
1071 int vmin = *( std::min_element( gids.begin(), gids.end() ) ),
1072 vmax = *( std::max_element( gids.begin(), gids.end() ) );
1073 dbgOut.tprintf( 1, "Vertex gids %d-%d\n", vmin, vmax );
1074 #endif
1075
1076 // Add elements to the range passed in
1077 faces.insert( scd_box->start_element(), scd_box->start_element() + scd_box->num_elements() - 1 );
1078
1079 if( 2 <= dbgOut.get_verbosity() )
1080 {
1081 assert( scd_box->boundary_complete() );
1082 EntityHandle dum_ent = scd_box->start_element();
1083 rval = mbImpl->list_entities( &dum_ent, 1 );MB_CHK_SET_ERR( rval, "Trouble listing first hex" );
1084
1085 std::vector< EntityHandle > connect;
1086 rval = mbImpl->get_connectivity( &dum_ent, 1, connect );MB_CHK_SET_ERR( rval, "Trouble getting connectivity" );
1087
1088 rval = mbImpl->list_entities( &connect[0], connect.size() );MB_CHK_SET_ERR( rval, "Trouble listing element connectivity" );
1089 }
1090
1091 Range edges;
1092 mbImpl->get_adjacencies( faces, 1, true, edges, Interface::UNION );
1093
1094 // Create COORDS tag for quads
1095 rval = create_quad_coordinate_tag();MB_CHK_SET_ERR( rval, "Trouble creating COORDS tag for quads" );
1096
1097 return MB_SUCCESS;
1098 }
1099
1100 ErrorCode ScdNCHelper::read_variables( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
1101 {
1102 std::vector< ReadNC::VarData > vdatas;
1103 std::vector< ReadNC::VarData > vsetdatas;
1104
1105 ErrorCode rval = read_variables_setup( var_names, tstep_nums, vdatas, vsetdatas );MB_CHK_SET_ERR( rval, "Trouble setting up to read variables" );
1106
1107 if( !vsetdatas.empty() )
1108 {
1109 rval = read_variables_to_set( vsetdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to set" );
1110 }
1111
1112 if( !vdatas.empty() )
1113 {
1114 rval = read_scd_variables_to_nonset( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to verts/edges/faces" );
1115 }
1116
1117 return MB_SUCCESS;
1118 }
1119
1120 ErrorCode ScdNCHelper::read_scd_variables_to_nonset_allocate( std::vector< ReadNC::VarData >& vdatas,
1121 std::vector< int >& tstep_nums )
1122 {
1123 Interface*& mbImpl = _readNC->mbImpl;
1124 std::vector< int >& dimLens = _readNC->dimLens;
1125 DebugOutput& dbgOut = _readNC->dbgOut;
1126
1127 Range* range = NULL;
1128
1129 // Get vertices
1130 Range verts;
1131 ErrorCode rval = mbImpl->get_entities_by_dimension( _fileSet, 0, verts );MB_CHK_SET_ERR( rval, "Trouble getting vertices in current file set" );
1132 assert( "Should only have a single vertex subrange, since they were read in one shot" && verts.psize() == 1 );
1133
1134 Range edges;
1135 rval = mbImpl->get_entities_by_dimension( _fileSet, 1, edges );MB_CHK_SET_ERR( rval, "Trouble getting edges in current file set" );
1136
1137 // Get faces
1138 Range faces;
1139 rval = mbImpl->get_entities_by_dimension( _fileSet, 2, faces );MB_CHK_SET_ERR( rval, "Trouble getting faces in current file set" );
1140 assert( "Should only have a single face subrange, since they were read in one shot" && faces.psize() == 1 );
1141
1142 #ifdef MOAB_HAVE_MPI
1143 moab::Range faces_owned;
1144 bool& isParallel = _readNC->isParallel;
1145 if( isParallel )
1146 {
1147 ParallelComm*& myPcomm = _readNC->myPcomm;
1148 rval = myPcomm->filter_pstatus( faces, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &faces_owned );MB_CHK_SET_ERR( rval, "Trouble getting owned faces in current file set" );
1149 }
1150 else
1151 faces_owned = faces; // Not running in parallel, but still with MPI
1152 #endif
1153
1154 for( unsigned int i = 0; i < vdatas.size(); i++ )
1155 {
1156 // Support non-set variables with 4 dimensions like (time, lev, lat, lon)
1157 assert( 4 == vdatas[i].varDims.size() );
1158
1159 // For a non-set variable, time should be the first dimension
1160 assert( tDim == vdatas[i].varDims[0] );
1161
1162 // Set up readStarts and readCounts
1163 vdatas[i].readStarts.resize( 4 );
1164 vdatas[i].readCounts.resize( 4 );
1165
1166 // First: time
1167 vdatas[i].readStarts[0] = 0; // This value is timestep dependent, will be set later
1168 vdatas[i].readCounts[0] = 1;
1169
1170 // Next: lev
1171 vdatas[i].readStarts[1] = 0;
1172 vdatas[i].readCounts[1] = vdatas[i].numLev;
1173
1174 // Finally: lat (or slat) and lon (or slon)
1175 switch( vdatas[i].entLoc )
1176 {
1177 case ReadNC::ENTLOCVERT:
1178 // Vertices
1179 vdatas[i].readStarts[2] = lDims[1];
1180 vdatas[i].readCounts[2] = lDims[4] - lDims[1] + 1;
1181 vdatas[i].readStarts[3] = lDims[0];
1182 vdatas[i].readCounts[3] = lDims[3] - lDims[0] + 1;
1183 range = &verts;
1184 break;
1185 case ReadNC::ENTLOCNSEDGE:
1186 case ReadNC::ENTLOCEWEDGE:
1187 case ReadNC::ENTLOCEDGE:
1188 // Not implemented yet, set a global error
1189 MB_SET_GLB_ERR( MB_NOT_IMPLEMENTED, "Reading edge data is not implemented yet" );
1190 case ReadNC::ENTLOCFACE:
1191 // Faces
1192 vdatas[i].readStarts[2] = lCDims[1];
1193 vdatas[i].readCounts[2] = lCDims[4] - lCDims[1] + 1;
1194 vdatas[i].readStarts[3] = lCDims[0];
1195 vdatas[i].readCounts[3] = lCDims[3] - lCDims[0] + 1;
1196 #ifdef MOAB_HAVE_MPI
1197 range = &faces_owned;
1198 #else
1199 range = &faces;
1200 #endif
1201 break;
1202 default:
1203 MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << vdatas[i].varName );
1204 }
1205
1206 for( unsigned int t = 0; t < tstep_nums.size(); t++ )
1207 {
1208 dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
1209
1210 if( tstep_nums[t] >= dimLens[tDim] )
1211 {
1212 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] );
1213 }
1214
1215 // Get the tag to read into
1216 if( !vdatas[i].varTags[t] )
1217 {
1218 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 to non-set variable " << vdatas[i].varName );
1219 }
1220
1221 // Get ptr to tag space
1222 void* data;
1223 int count;
1224 rval = mbImpl->tag_iterate( vdatas[i].varTags[t], range->begin(), range->end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate tag for non-set variable " << vdatas[i].varName );
1225 assert( (unsigned)count == range->size() );
1226 vdatas[i].varDatas[t] = data;
1227 }
1228
1229 // Get variable size
1230 vdatas[i].sz = 1;
1231 for( std::size_t idx = 0; idx != vdatas[i].readCounts.size(); idx++ )
1232 vdatas[i].sz *= vdatas[i].readCounts[idx];
1233 }
1234
1235 return rval;
1236 }
1237
1238 ErrorCode ScdNCHelper::read_scd_variables_to_nonset( std::vector< ReadNC::VarData >& vdatas,
1239 std::vector< int >& tstep_nums )
1240 {
1241 DebugOutput& dbgOut = _readNC->dbgOut;
1242
1243 ErrorCode rval = read_scd_variables_to_nonset_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read non-set variables" );
1244
1245 // Finally, read into that space
1246 int success;
1247 for( unsigned int i = 0; i < vdatas.size(); i++ )
1248 {
1249 std::size_t sz = vdatas[i].sz;
1250
1251 // A typical supported variable: float T(time, lev, lat, lon)
1252 // For tag values, need transpose (lev, lat, lon) to (lat, lon, lev)
1253 size_t ni = vdatas[i].readCounts[3]; // lon or slon
1254 size_t nj = vdatas[i].readCounts[2]; // lat or slat
1255 size_t nk = vdatas[i].readCounts[1]; // lev
1256
1257 for( unsigned int t = 0; t < tstep_nums.size(); t++ )
1258 {
1259 // Tag data for this timestep
1260 void* data = vdatas[i].varDatas[t];
1261
1262 // Set readStart for each timestep along time dimension
1263 vdatas[i].readStarts[0] = tstep_nums[t];
1264
1265 switch( vdatas[i].varDataType )
1266 {
1267 case NC_BYTE:
1268 case NC_CHAR: {
1269 std::vector< char > tmpchardata( sz );
1270 success = NCFUNCAG( _vara_text )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
1271 &vdatas[i].readCounts[0], &tmpchardata[0] );
1272 if( success )
1273 MB_SET_ERR( MB_FAILURE, "Failed to read byte/char data for variable " << vdatas[i].varName );
1274 if( vdatas[i].numLev > 1 )
1275 // Transpose (lev, lat, lon) to (lat, lon, lev)
1276 kji_to_jik( ni, nj, nk, data, &tmpchardata[0] );
1277 else
1278 {
1279 for( std::size_t idx = 0; idx != tmpchardata.size(); idx++ )
1280 ( (char*)data )[idx] = tmpchardata[idx];
1281 }
1282 break;
1283 }
1284 case NC_SHORT:
1285 case NC_INT: {
1286 std::vector< int > tmpintdata( sz );
1287 success = NCFUNCAG( _vara_int )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
1288 &vdatas[i].readCounts[0], &tmpintdata[0] );
1289 if( success )
1290 MB_SET_ERR( MB_FAILURE, "Failed to read short/int data for variable " << vdatas[i].varName );
1291 if( vdatas[i].numLev > 1 )
1292 // Transpose (lev, lat, lon) to (lat, lon, lev)
1293 kji_to_jik( ni, nj, nk, data, &tmpintdata[0] );
1294 else
1295 {
1296 for( std::size_t idx = 0; idx != tmpintdata.size(); idx++ )
1297 ( (int*)data )[idx] = tmpintdata[idx];
1298 }
1299 break;
1300 }
1301 case NC_INT64: {
1302 std::vector< long > tmpintdata( sz );
1303 success = NCFUNCAG( _vara_long )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
1304 &vdatas[i].readCounts[0], &tmpintdata[0] );
1305 if( success )
1306 MB_SET_ERR( MB_FAILURE, "Failed to read long data for variable " << vdatas[i].varName );
1307 if( vdatas[i].numLev > 1 )
1308 // Transpose (lev, lat, lon) to (lat, lon, lev)
1309 kji_to_jik( ni, nj, nk, data, &tmpintdata[0] );
1310 else
1311 {
1312 for( std::size_t idx = 0; idx != tmpintdata.size(); idx++ )
1313 ( (long*)data )[idx] = tmpintdata[idx];
1314 }
1315 break;
1316 }
1317 case NC_FLOAT:
1318 case NC_DOUBLE: {
1319 std::vector< double > tmpdoubledata( sz );
1320 success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
1321 &vdatas[i].readCounts[0], &tmpdoubledata[0] );
1322 if( success )
1323 MB_SET_ERR( MB_FAILURE, "Failed to read float/double data for variable " << vdatas[i].varName );
1324 if( vdatas[i].numLev > 1 )
1325 // Transpose (lev, lat, lon) to (lat, lon, lev)
1326 kji_to_jik( ni, nj, nk, data, &tmpdoubledata[0] );
1327 else
1328 {
1329 for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
1330 ( (double*)data )[idx] = tmpdoubledata[idx];
1331 }
1332 break;
1333 }
1334 default:
1335 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
1336 }
1337 }
1338 }
1339
1340 // Debug output, if requested
1341 if( 1 == dbgOut.get_verbosity() )
1342 {
1343 dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
1344 for( unsigned int i = 1; i < vdatas.size(); i++ )
1345 dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
1346 dbgOut.tprintf( 1, "\n" );
1347 }
1348
1349 return rval;
1350 }
1351
1352 ErrorCode ScdNCHelper::create_quad_coordinate_tag()
1353 {
1354 Interface*& mbImpl = _readNC->mbImpl;
1355
1356 Range ents;
1357 ErrorCode rval = mbImpl->get_entities_by_type( _fileSet, moab::MBQUAD, ents );MB_CHK_SET_ERR( rval, "Trouble getting quads" );
1358
1359 std::size_t numOwnedEnts = 0;
1360 #ifdef MOAB_HAVE_MPI
1361 Range ents_owned;
1362 bool& isParallel = _readNC->isParallel;
1363 if( isParallel )
1364 {
1365 ParallelComm*& myPcomm = _readNC->myPcomm;
1366 rval = myPcomm->filter_pstatus( ents, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &ents_owned );MB_CHK_SET_ERR( rval, "Trouble getting owned quads" );
1367 numOwnedEnts = ents_owned.size();
1368 }
1369 else
1370 {
1371 numOwnedEnts = ents.size();
1372 ents_owned = ents;
1373 }
1374 #else
1375 numOwnedEnts = ents.size();
1376 #endif
1377
1378 if( numOwnedEnts == 0 ) return MB_SUCCESS;
1379
1380 assert( numOwnedEnts == ilCVals.size() * jlCVals.size() );
1381 std::vector< double > coords( numOwnedEnts * 3 );
1382 std::size_t pos = 0;
1383 for( std::size_t j = 0; j != jlCVals.size(); ++j )
1384 {
1385 for( std::size_t i = 0; i != ilCVals.size(); ++i )
1386 {
1387 pos = j * ilCVals.size() * 3 + i * 3;
1388 coords[pos] = ilCVals[i];
1389 coords[pos + 1] = jlCVals[j];
1390 coords[pos + 2] = 0.0;
1391 }
1392 }
1393 std::string tag_name = "COORDS";
1394 Tag tagh = 0;
1395 rval = mbImpl->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating COORDS tag" );
1396
1397 void* data;
1398 int count;
1399 #ifdef MOAB_HAVE_MPI
1400 rval = mbImpl->tag_iterate( tagh, ents_owned.begin(), ents_owned.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate COORDS tag on quads" );
1401 #else
1402 rval = mbImpl->tag_iterate( tagh, ents.begin(), ents.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate COORDS tag on quads" );
1403 #endif
1404 assert( count == (int)numOwnedEnts );
1405 double* quad_data = (double*)data;
1406 std::copy( coords.begin(), coords.end(), quad_data );
1407
1408 return MB_SUCCESS;
1409 }
1410
1411 ErrorCode UcdNCHelper::read_variables( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
1412 {
1413 std::vector< ReadNC::VarData > vdatas;
1414 std::vector< ReadNC::VarData > vsetdatas;
1415
1416 ErrorCode rval = read_variables_setup( var_names, tstep_nums, vdatas, vsetdatas );MB_CHK_SET_ERR( rval, "Trouble setting up to read variables" );
1417
1418 if( !vsetdatas.empty() )
1419 {
1420 rval = read_variables_to_set( vsetdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to set" );
1421 }
1422
1423 if( !vdatas.empty() )
1424 {
1425 #ifdef MOAB_HAVE_PNETCDF
1426 // With pnetcdf support, we will use async read
1427 rval = read_ucd_variables_to_nonset_async( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to verts/edges/faces" );
1428 #else
1429 // Without pnetcdf support, we will use old read
1430 rval = read_ucd_variables_to_nonset( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to verts/edges/faces" );
1431 #endif
1432 }
1433
1434 return MB_SUCCESS;
1435 }
1436
1437 static double tolerance = 1.e-12;
1438 // Used to put points in an STL tree-based container
1439 bool NCHelper::Node3D::operator<( const NCHelper::Node3D& other ) const
1440 {
1441 if( coords[0] <= other.coords[0] - tolerance )
1442 return true;
1443 else if( coords[0] >= other.coords[0] + tolerance )
1444 return false;
1445 if( coords[1] <= other.coords[1] - tolerance )
1446 return true;
1447 else if( coords[1] >= other.coords[1] + tolerance )
1448 return false;
1449 if( coords[2] <= other.coords[2] - tolerance )
1450 return true;
1451 else if( coords[0] >= other.coords[0] + tolerance )
1452 return false;
1453
1454 return false;
1455 }
1456 /*bool NCHelper::Node3D::operator == ( const NCHelper::Node3D& other ) const
1457 return ( !(*this < other) && ! ( other < *this) ) ;*/
1458 } // namespace moab
1.9.1.