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