 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
00001 #include "NCHelper.hpp"
00002 #include "NCHelperEuler.hpp"
00003 #include "NCHelperFV.hpp"
00004 #include "NCHelperDomain.hpp"
00005 #include "NCHelperScrip.hpp"
00006 #include "NCHelperHOMME.hpp"
00007 #include "NCHelperMPAS.hpp"
00008 #include "NCHelperGCRM.hpp"
00009
00010 #include
00011
00012 #include "MBTagConventions.hpp"
00013
00014 #ifdef WIN32
00015 #ifdef size_t
00016 #undef size_t
00017 #endif
00018 #endif
00019
00020 namespace moab
00021 {
00022
00023 NCHelper* NCHelper::get_nc_helper( ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet )
00024 {
00025 // Check if CF convention is being followed
00026 bool is_CF = false;
00027
00028 std::map< std::string, ReadNC::AttData >& globalAtts = readNC->globalAtts;
00029 std::map< std::string, ReadNC::AttData >::iterator attIt = globalAtts.find( "conventions" );
00030 if( attIt == globalAtts.end() ) attIt = globalAtts.find( "Conventions" );
00031
00032 if( attIt != globalAtts.end() )
00033 {
00034 unsigned int sz = attIt->second.attLen;
00035 std::string att_data;
00036 att_data.resize( sz + 1 );
00037 att_data[sz] = '\000';
00038 int success =
00039 NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
00040 if( 0 == success && att_data.find( "CF" ) != std::string::npos ) is_CF = true;
00041 }
00042
00043 if( is_CF )
00044 {
00045 if( NCHelperEuler::can_read_file( readNC, fileId ) )
00046 return new( std::nothrow ) NCHelperEuler( readNC, fileId, opts, fileSet );
00047 else if( NCHelperFV::can_read_file( readNC, fileId ) )
00048 return new( std::nothrow ) NCHelperFV( readNC, fileId, opts, fileSet );
00049 else if( NCHelperHOMME::can_read_file( readNC, fileId ) )
00050 return new( std::nothrow ) NCHelperHOMME( readNC, fileId, opts, fileSet );
00051 else if( NCHelperDomain::can_read_file( readNC, fileId ) )
00052 return new( std::nothrow ) NCHelperDomain( readNC, fileId, opts, fileSet );
00053 }
00054 else
00055 {
00056 if( NCHelperMPAS::can_read_file( readNC ) )
00057 return new( std::nothrow ) NCHelperMPAS( readNC, fileId, opts, fileSet );
00058 // For a HOMME connectivity file, there might be no CF convention
00059 else if( NCHelperHOMME::can_read_file( readNC, fileId ) )
00060 return new( std::nothrow ) NCHelperHOMME( readNC, fileId, opts, fileSet );
00061 // gcrm reader
00062 else if( NCHelperGCRM::can_read_file( readNC ) )
00063 return new( std::nothrow ) NCHelperGCRM( readNC, fileId, opts, fileSet );
00064 }
00065 // SCRIP can be CF or non CF, if it comes from MPAS :)
00066 if( NCHelperScrip::can_read_file( readNC, fileId ) )
00067 return new( std::nothrow ) NCHelperScrip( readNC, fileId, opts, fileSet );
00068
00069 // Unknown NetCDF grid (will fill this in later for POP, CICE and CLM)
00070 return NULL;
00071 }
00072
00073 ErrorCode NCHelper::create_conventional_tags( const std::vector< int >& tstep_nums )
00074 {
00075 Interface*& mbImpl = _readNC->mbImpl;
00076 std::vector< std::string >& dimNames = _readNC->dimNames;
00077 std::vector< int >& dimLens = _readNC->dimLens;
00078 std::map< std::string, ReadNC::AttData >& globalAtts = _readNC->globalAtts;
00079 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
00080 DebugOutput& dbgOut = _readNC->dbgOut;
00081 int& partMethod = _readNC->partMethod;
00082 ScdInterface* scdi = _readNC->scdi;
00083
00084 ErrorCode rval;
00085 std::string tag_name;
00086
00087 // <__NUM_DIMS>
00088 Tag numDimsTag = 0;
00089 tag_name = "__NUM_DIMS";
00090 int numDims = dimNames.size();
00091 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 );
00092 rval = mbImpl->tag_set_data( numDimsTag, &_fileSet, 1, &numDims );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00093 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00094
00095 // <__NUM_VARS>
00096 Tag numVarsTag = 0;
00097 tag_name = "__NUM_VARS";
00098 int numVars = varInfo.size();
00099 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 );
00100 rval = mbImpl->tag_set_data( numVarsTag, &_fileSet, 1, &numVars );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00101 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00102
00103 // <__DIM_NAMES>
00104 Tag dimNamesTag = 0;
00105 tag_name = "__DIM_NAMES";
00106 std::string dimnames;
00107 unsigned int dimNamesSz = dimNames.size();
00108 for( unsigned int i = 0; i != dimNamesSz; i++ )
00109 {
00110 dimnames.append( dimNames[i] );
00111 dimnames.push_back( '\0' );
00112 }
00113 int dimnamesSz = dimnames.size();
00114 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, dimNamesTag,
00115 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00116 const void* ptr = dimnames.c_str();
00117 rval = mbImpl->tag_set_by_ptr( dimNamesTag, &_fileSet, 1, &ptr, &dimnamesSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00118 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00119
00120 // <__DIM_LENS>
00121 Tag dimLensTag = 0;
00122 tag_name = "__DIM_LENS";
00123 int dimLensSz = dimLens.size();
00124 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, dimLensTag,
00125 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00126 ptr = &( dimLens[0] );
00127 rval = mbImpl->tag_set_by_ptr( dimLensTag, &_fileSet, 1, &ptr, &dimLensSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00128 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00129
00130 // <__VAR_NAMES>
00131 Tag varNamesTag = 0;
00132 tag_name = "__VAR_NAMES";
00133 std::string varnames;
00134 std::map< std::string, ReadNC::VarData >::iterator mapIter;
00135 for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
00136 {
00137 varnames.append( mapIter->first );
00138 varnames.push_back( '\0' );
00139 }
00140 int varnamesSz = varnames.size();
00141 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varNamesTag,
00142 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00143 ptr = varnames.c_str();
00144 rval = mbImpl->tag_set_by_ptr( varNamesTag, &_fileSet, 1, &ptr, &varnamesSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00145 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00146
00147 // ___LOC_MINMAX (for time)
00148 for( unsigned int i = 0; i != dimNamesSz; i++ )
00149 {
00150 if( dimNames[i] == "time" || dimNames[i] == "Time" || dimNames[i] == "t" )
00151 {
00152 // some files might have Time dimension as 0, it will not appear in any
00153 // variables, so skip it
00154 if( nTimeSteps == 0 ) continue;
00155 std::stringstream ss_tag_name;
00156 ss_tag_name << "__" << dimNames[i] << "_LOC_MINMAX";
00157 tag_name = ss_tag_name.str();
00158 Tag tagh = 0;
00159 std::vector< int > val( 2, 0 );
00160 val[0] = 0;
00161 val[1] = nTimeSteps - 1;
00162 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 );
00163 rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00164 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00165 }
00166 }
00167
00168 // ___LOC_VALS (for time)
00169 for( unsigned int i = 0; i != dimNamesSz; i++ )
00170 {
00171 if( dimNames[i] == "time" || dimNames[i] == "Time" || dimNames[i] == "t" )
00172 {
00173 std::vector< int > val;
00174 if( !tstep_nums.empty() )
00175 val = tstep_nums;
00176 else
00177 {
00178 // some files might have Time dimension as 0, it will not appear in any
00179 // variables, so skip it
00180 if( tVals.empty() ) continue;
00181 val.resize( tVals.size() );
00182 for( unsigned int j = 0; j != tVals.size(); j++ )
00183 val[j] = j;
00184 }
00185 Tag tagh = 0;
00186 std::stringstream ss_tag_name;
00187 ss_tag_name << "__" << dimNames[i] << "_LOC_VALS";
00188 tag_name = ss_tag_name.str();
00189 rval = mbImpl->tag_get_handle( tag_name.c_str(), val.size(), MB_TYPE_INTEGER, tagh,
00190 MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00191 rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00192 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00193 }
00194 }
00195
00196 // ___DIMS
00197 for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
00198 {
00199 Tag varNamesDimsTag = 0;
00200 std::stringstream ss_tag_name;
00201 ss_tag_name << "__" << mapIter->first << "_DIMS";
00202 tag_name = ss_tag_name.str();
00203 unsigned int varDimSz = varInfo[mapIter->first].varDims.size();
00204 if( varDimSz == 0 ) continue;
00205 std::vector< Tag > varDimTags( varDimSz );
00206 for( unsigned int i = 0; i != varDimSz; i++ )
00207 {
00208 Tag tmptag = 0;
00209 std::string tmptagname = dimNames[varInfo[mapIter->first].varDims[i]];
00210 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 );
00211 varDimTags[i] = tmptag;
00212 }
00213 // rval = mbImpl->tag_get_handle(tag_name.c_str(), varDimSz, MB_TYPE_HANDLE,
00214 // varNamesDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT); We should not use MB_TYPE_HANDLE for Tag
00215 // here. Tag is a pointer, which is 4 bytes on 32 bit machines and 8 bytes on 64 bit
00216 // machines. Normally, entity handle is 8 bytes on 64 bit machines, but it can also be
00217 // configured to 4 bytes.
00218 rval = mbImpl->tag_get_handle( tag_name.c_str(), varDimSz * sizeof( Tag ), MB_TYPE_OPAQUE, varNamesDimsTag,
00219 MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00220 rval = mbImpl->tag_set_data( varNamesDimsTag, &_fileSet, 1, &( varDimTags[0] ) );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00221 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00222 }
00223
00224 //
00225 Tag part_tag = scdi->part_method_tag();
00226 if( !part_tag ) MB_SET_ERR( MB_FAILURE, "Trouble getting PARTITION_METHOD tag" );
00227 rval = mbImpl->tag_set_data( part_tag, &_fileSet, 1, &partMethod );MB_CHK_SET_ERR( rval, "Trouble setting data to PARTITION_METHOD tag" );
00228 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00229
00230 // <__GLOBAL_ATTRIBS>
00231 tag_name = "__GLOBAL_ATTRIBS";
00232 Tag globalAttTag = 0;
00233 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, globalAttTag,
00234 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00235 std::string gattVal;
00236 std::vector< int > gattLen;
00237 rval = create_attrib_string( globalAtts, gattVal, gattLen );MB_CHK_SET_ERR( rval, "Trouble creating global attribute string" );
00238 const void* gattptr = gattVal.c_str();
00239 int globalAttSz = gattVal.size();
00240 rval = mbImpl->tag_set_by_ptr( globalAttTag, &_fileSet, 1, &gattptr, &globalAttSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00241 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00242
00243 // <__GLOBAL_ATTRIBS_LEN>
00244 tag_name = "__GLOBAL_ATTRIBS_LEN";
00245 Tag globalAttLenTag = 0;
00246 if( gattLen.size() == 0 ) gattLen.push_back( 0 );
00247 rval = mbImpl->tag_get_handle( tag_name.c_str(), gattLen.size(), MB_TYPE_INTEGER, globalAttLenTag,
00248 MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00249 rval = mbImpl->tag_set_data( globalAttLenTag, &_fileSet, 1, &gattLen[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00250 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00251
00252 // ___ATTRIBS and ___ATTRIBS_LEN
00253 for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
00254 {
00255 std::stringstream ssTagName;
00256 ssTagName << "__" << mapIter->first << "_ATTRIBS";
00257 tag_name = ssTagName.str();
00258 Tag varAttTag = 0;
00259 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varAttTag,
00260 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00261
00262 std::string varAttVal;
00263 std::vector< int > varAttLen;
00264 if( mapIter->second.numAtts < 1 )
00265 {
00266 if( dummyVarNames.find( mapIter->first ) != dummyVarNames.end() )
00267 {
00268 // This variable is a dummy coordinate variable
00269 varAttVal = "DUMMY_VAR";
00270 }
00271 else
00272 {
00273 // This variable has no attributes
00274 varAttVal = "NO_ATTRIBS";
00275 }
00276 }
00277 else
00278 {
00279 rval = create_attrib_string( mapIter->second.varAtts, varAttVal, varAttLen );MB_CHK_SET_ERR( rval, "Trouble creating attribute string for variable " << mapIter->first );
00280 }
00281 const void* varAttPtr = varAttVal.c_str();
00282 int varAttSz = varAttVal.size();
00283 if( 0 == varAttSz ) varAttSz = 1;
00284 rval = mbImpl->tag_set_by_ptr( varAttTag, &_fileSet, 1, &varAttPtr, &varAttSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00285 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00286
00287 ssTagName << "_LEN";
00288 tag_name = ssTagName.str();
00289 Tag varAttLenTag = 0;
00290 if( 0 == varAttLen.size() ) varAttLen.push_back( 0 );
00291 rval = mbImpl->tag_get_handle( tag_name.c_str(), varAttLen.size(), MB_TYPE_INTEGER, varAttLenTag,
00292 MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00293 rval = mbImpl->tag_set_data( varAttLenTag, &_fileSet, 1, &varAttLen[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00294 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00295 }
00296
00297 // <__VAR_NAMES_LOCATIONS>
00298 tag_name = "__VAR_NAMES_LOCATIONS";
00299 Tag varNamesLocsTag = 0;
00300 std::vector< int > varNamesLocs( varInfo.size() );
00301 rval = mbImpl->tag_get_handle( tag_name.c_str(), varNamesLocs.size(), MB_TYPE_INTEGER, varNamesLocsTag,
00302 MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00303 for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
00304 {
00305 varNamesLocs[std::distance( varInfo.begin(), mapIter )] = mapIter->second.entLoc;
00306 }
00307 rval = mbImpl->tag_set_data( varNamesLocsTag, &_fileSet, 1, &varNamesLocs[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00308 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00309
00310 // <__MESH_TYPE>
00311 Tag meshTypeTag = 0;
00312 tag_name = "__MESH_TYPE";
00313 std::string meshTypeName = get_mesh_type_name();
00314
00315 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, meshTypeTag,
00316 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00317 ptr = meshTypeName.c_str();
00318 int leng = meshTypeName.size();
00319 rval = mbImpl->tag_set_by_ptr( meshTypeTag, &_fileSet, 1, &ptr, &leng );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00320 dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
00321
00322 return MB_SUCCESS;
00323 }
00324
00325 ErrorCode NCHelper::update_time_tag_vals()
00326 {
00327 Interface*& mbImpl = _readNC->mbImpl;
00328 std::vector< std::string >& dimNames = _readNC->dimNames;
00329
00330 ErrorCode rval;
00331
00332 // The time tag might be a dummy one (e.g. 'Time' for MPAS)
00333 std::string time_tag_name = dimNames[tDim];
00334 if( dummyVarNames.find( time_tag_name ) != dummyVarNames.end() ) return MB_SUCCESS;
00335
00336 Tag time_tag = 0;
00337 const void* data = NULL;
00338 int time_tag_size = 0;
00339 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 );
00340 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 );
00341 const double* time_tag_vals = static_cast< const double* >( data );
00342
00343 // Merge tVals (read from current file) to existing time tag
00344 // Assume that time_tag_vals and tVals are both sorted
00345 std::vector< double > merged_time_vals;
00346 merged_time_vals.reserve( time_tag_size + nTimeSteps );
00347 int i = 0;
00348 int j = 0;
00349
00350 // Merge time values from time_tag_vals and tVals
00351 while( i < time_tag_size && j < nTimeSteps )
00352 {
00353 if( time_tag_vals[i] < tVals[j] )
00354 merged_time_vals.push_back( time_tag_vals[i++] );
00355 else
00356 merged_time_vals.push_back( tVals[j++] );
00357 }
00358
00359 // Append remaining time values of time_tag_vals (if any)
00360 while( i < time_tag_size )
00361 merged_time_vals.push_back( time_tag_vals[i++] );
00362
00363 // Append remaining time values of tVals (if any)
00364 while( j < nTimeSteps )
00365 merged_time_vals.push_back( tVals[j++] );
00366
00367 data = &merged_time_vals[0];
00368 time_tag_size = merged_time_vals.size();
00369 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 );
00370
00371 return MB_SUCCESS;
00372 }
00373
00374 ErrorCode NCHelper::read_variables_setup( std::vector< std::string >& var_names,
00375 std::vector< int >& tstep_nums,
00376 std::vector< ReadNC::VarData >& vdatas,
00377 std::vector< ReadNC::VarData >& vsetdatas )
00378 {
00379 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
00380 std::vector< std::string >& dimNames = _readNC->dimNames;
00381
00382 std::map< std::string, ReadNC::VarData >::iterator mit;
00383
00384 // If empty read them all (except ignored variables)
00385 if( var_names.empty() )
00386 {
00387 for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
00388 {
00389 ReadNC::VarData vd = ( *mit ).second;
00390
00391 // If read all variables at once, skip ignored ones
00392 if( ignoredVarNames.find( vd.varName ) != ignoredVarNames.end() ) continue;
00393
00394 // Coordinate variables (include dummy ones) were read to the file set by default
00395 if( std::find( dimNames.begin(), dimNames.end(), vd.varName ) != dimNames.end() ) continue;
00396
00397 if( vd.entLoc == ReadNC::ENTLOCSET )
00398 vsetdatas.push_back( vd );
00399 else
00400 vdatas.push_back( vd );
00401 }
00402 }
00403 else
00404 {
00405 // Read specified variables (might include ignored ones)
00406 for( unsigned int i = 0; i < var_names.size(); i++ )
00407 {
00408 mit = varInfo.find( var_names[i] );
00409 if( mit != varInfo.end() )
00410 {
00411 ReadNC::VarData vd = ( *mit ).second;
00412
00413 // Upon creation of dummy coordinate variables, tag values have already been set
00414 if( dummyVarNames.find( vd.varName ) != dummyVarNames.end() ) continue;
00415
00416 if( vd.entLoc == ReadNC::ENTLOCSET )
00417 vsetdatas.push_back( vd );
00418 else
00419 vdatas.push_back( vd );
00420 }
00421 else
00422 {
00423 MB_SET_ERR( MB_FAILURE, "Couldn't find specified variable " << var_names[i] );
00424 }
00425 }
00426 }
00427
00428 if( tstep_nums.empty() && nTimeSteps > 0 )
00429 {
00430 // No timesteps input, get them all
00431 for( int i = 0; i < nTimeSteps; i++ )
00432 tstep_nums.push_back( i );
00433 }
00434
00435 if( !tstep_nums.empty() )
00436 {
00437 for( unsigned int i = 0; i < vdatas.size(); i++ )
00438 {
00439 vdatas[i].varTags.resize( tstep_nums.size(), 0 );
00440 vdatas[i].varDatas.resize( tstep_nums.size() );
00441 // NC reader assumes that non-set variables always have timesteps
00442 assert( std::find( vdatas[i].varDims.begin(), vdatas[i].varDims.end(), tDim ) != vdatas[i].varDims.end() );
00443 vdatas[i].has_tsteps = true;
00444 }
00445
00446 for( unsigned int i = 0; i < vsetdatas.size(); i++ )
00447 {
00448 if( ( std::find( vsetdatas[i].varDims.begin(), vsetdatas[i].varDims.end(), tDim ) !=
00449 vsetdatas[i].varDims.end() ) &&
00450 ( vsetdatas[i].varName != dimNames[tDim] ) )
00451 {
00452 // Set variables with timesteps: e.g. xtime(Time) or xtime(Time, StrLen)
00453 vsetdatas[i].varTags.resize( tstep_nums.size(), 0 );
00454 vsetdatas[i].varDatas.resize( tstep_nums.size() );
00455 vsetdatas[i].has_tsteps = true;
00456 }
00457 else
00458 {
00459 // Set variables without timesteps: no time dimension, or time itself
00460 vsetdatas[i].varTags.resize( 1, 0 );
00461 vsetdatas[i].varDatas.resize( 1 );
00462 vsetdatas[i].has_tsteps = false;
00463 }
00464 }
00465 }
00466
00467 return MB_SUCCESS;
00468 }
00469
00470 ErrorCode NCHelper::read_variables_to_set( std::vector< ReadNC::VarData >& vdatas, std::vector< int >& tstep_nums )
00471 {
00472 Interface*& mbImpl = _readNC->mbImpl;
00473 DebugOutput& dbgOut = _readNC->dbgOut;
00474
00475 ErrorCode rval = read_variables_to_set_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read set variables" );
00476
00477 // Finally, read into that space
00478 int success;
00479 for( unsigned int i = 0; i < vdatas.size(); i++ )
00480 {
00481 // Note, for set variables without timesteps, loop one time and then break
00482 for( unsigned int t = 0; t < tstep_nums.size(); t++ )
00483 {
00484 void* data = vdatas[i].varDatas[t];
00485
00486 // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen)
00487 if( vdatas[i].has_tsteps )
00488 {
00489 // Set readStart for each timestep along time dimension
00490 vdatas[i].readStarts[0] = tstep_nums[t];
00491 }
00492
00493 switch( vdatas[i].varDataType )
00494 {
00495 case NC_BYTE:
00496 case NC_CHAR:
00497 success = NCFUNCAG( _vara_text )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
00498 &vdatas[i].readCounts[0], (char*)data );
00499 if( success )
00500 MB_SET_ERR( MB_FAILURE, "Failed to read byte/char data for variable " << vdatas[i].varName );
00501 break;
00502 case NC_SHORT:
00503 case NC_INT:
00504 success = NCFUNCAG( _vara_int )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
00505 &vdatas[i].readCounts[0], (int*)data );
00506 if( success )
00507 MB_SET_ERR( MB_FAILURE, "Failed to read short/int data for variable " << vdatas[i].varName );
00508 break;
00509 case NC_FLOAT:
00510 case NC_DOUBLE:
00511 success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
00512 &vdatas[i].readCounts[0], (double*)data );
00513 if( success )
00514 MB_SET_ERR( MB_FAILURE, "Failed to read float/double data for variable " << vdatas[i].varName );
00515 break;
00516 default:
00517 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
00518 }
00519
00520 dbgOut.tprintf( 2, "Setting data for variable %s, time step %d\n", vdatas[i].varName.c_str(),
00521 tstep_nums[t] );
00522 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 );
00523
00524 // Memory pointed by pointer data can be deleted, as tag_set_by_ptr() has already copied
00525 // the tag values
00526 switch( vdatas[i].varDataType )
00527 {
00528 case NC_BYTE:
00529 case NC_CHAR:
00530 delete[]( char* ) data;
00531 break;
00532 case NC_SHORT:
00533 case NC_INT:
00534 delete[]( int* ) data;
00535 break;
00536 case NC_FLOAT:
00537 case NC_DOUBLE:
00538 delete[]( double* ) data;
00539 break;
00540 default:
00541 break;
00542 }
00543 vdatas[i].varDatas[t] = NULL;
00544
00545 // Loop continues only for set variables with timesteps, e.g. xtime(Time) or xtime(Time,
00546 // StrLen)
00547 if( !vdatas[i].has_tsteps ) break;
00548 }
00549 }
00550
00551 // Debug output, if requested
00552 if( 1 == dbgOut.get_verbosity() )
00553 {
00554 dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
00555 for( unsigned int i = 1; i < vdatas.size(); i++ )
00556 dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
00557 dbgOut.tprintf( 1, "\n" );
00558 }
00559
00560 return rval;
00561 }
00562
00563 ErrorCode NCHelper::read_coordinate( const char* var_name, int lmin, int lmax, std::vector< double >& cvals )
00564 {
00565 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
00566 std::map< std::string, ReadNC::VarData >::iterator vmit = varInfo.find( var_name );
00567 if( varInfo.end() == vmit ) MB_SET_ERR( MB_FAILURE, "Couldn't find variable " << var_name );
00568
00569 assert( lmin >= 0 && lmax >= lmin );
00570 NCDF_SIZE tstart = lmin;
00571 NCDF_SIZE tcount = lmax - lmin + 1;
00572 NCDF_DIFF dum_stride = 1;
00573 int success;
00574
00575 // Check size
00576 if( (std::size_t)tcount != cvals.size() ) cvals.resize( tcount );
00577
00578 // Check to make sure it's a float or double
00579 switch( ( *vmit ).second.varDataType )
00580 {
00581 case NC_FLOAT:
00582 case NC_DOUBLE:
00583 // Read float as double
00584 success =
00585 NCFUNCAG( _vars_double )( _fileId, ( *vmit ).second.varId, &tstart, &tcount, &dum_stride, &cvals[0] );
00586 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read float/double data for variable " << var_name );
00587 break;
00588 default:
00589 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << var_name );
00590 }
00591
00592 return MB_SUCCESS;
00593 }
00594
00595 ErrorCode NCHelper::get_tag_to_set( ReadNC::VarData& var_data, int tstep_num, Tag& tagh )
00596 {
00597 Interface*& mbImpl = _readNC->mbImpl;
00598 DebugOutput& dbgOut = _readNC->dbgOut;
00599 int& tStepBase = _readNC->tStepBase;
00600
00601 if( tStepBase > 0 ) tstep_num += tStepBase;
00602
00603 std::ostringstream tag_name;
00604 if( var_data.has_tsteps )
00605 tag_name << var_data.varName << tstep_num;
00606 else
00607 tag_name << var_data.varName;
00608
00609 ErrorCode rval = MB_SUCCESS;
00610 tagh = 0;
00611 switch( var_data.varDataType )
00612 {
00613 case NC_BYTE:
00614 case NC_CHAR:
00615 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), 0, MB_TYPE_OPAQUE, tagh,
00616 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
00617 break;
00618 case NC_SHORT:
00619 case NC_INT:
00620 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), 0, MB_TYPE_INTEGER, tagh,
00621 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
00622 break;
00623 case NC_FLOAT:
00624 case NC_DOUBLE:
00625 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), 0, MB_TYPE_DOUBLE, tagh,
00626 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
00627 break;
00628 default:
00629 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << var_data.varName );
00630 }
00631
00632 dbgOut.tprintf( 2, "Tag %s created\n", tag_name.str().c_str() );
00633
00634 return rval;
00635 }
00636
00637 ErrorCode NCHelper::get_tag_to_nonset( ReadNC::VarData& var_data, int tstep_num, Tag& tagh, int num_lev )
00638 {
00639 Interface*& mbImpl = _readNC->mbImpl;
00640 DebugOutput& dbgOut = _readNC->dbgOut;
00641 int& tStepBase = _readNC->tStepBase;
00642
00643 if( tStepBase > 0 ) tstep_num += tStepBase;
00644
00645 std::ostringstream tag_name;
00646 tag_name << var_data.varName << tstep_num;
00647
00648 ErrorCode rval = MB_SUCCESS;
00649 tagh = 0;
00650 switch( var_data.varDataType )
00651 {
00652 case NC_BYTE:
00653 case NC_CHAR:
00654 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), num_lev, MB_TYPE_OPAQUE, tagh,
00655 MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
00656 break;
00657 case NC_SHORT:
00658 case NC_INT:
00659 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), num_lev, MB_TYPE_INTEGER, tagh,
00660 MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
00661 break;
00662 case NC_FLOAT:
00663 case NC_DOUBLE:
00664 rval = mbImpl->tag_get_handle( tag_name.str().c_str(), num_lev, MB_TYPE_DOUBLE, tagh,
00665 MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
00666 break;
00667 default:
00668 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << var_data.varName );
00669 }
00670
00671 dbgOut.tprintf( 2, "Tag %s created\n", tag_name.str().c_str() );
00672
00673 return rval;
00674 }
00675
00676 ErrorCode NCHelper::create_attrib_string( const std::map< std::string, ReadNC::AttData >& attMap,
00677 std::string& attVal,
00678 std::vector< int >& attLen )
00679 {
00680 int success;
00681 std::stringstream ssAtt;
00682 unsigned int sz = 0;
00683 std::map< std::string, ReadNC::AttData >::const_iterator attIt = attMap.begin();
00684 for( ; attIt != attMap.end(); ++attIt )
00685 {
00686 ssAtt << attIt->second.attName;
00687 ssAtt << '\0';
00688 void* attData = NULL;
00689 switch( attIt->second.attDataType )
00690 {
00691 case NC_BYTE:
00692 case NC_CHAR:
00693 sz = attIt->second.attLen;
00694 attData = (char*)malloc( sz );
00695 success = NCFUNC( get_att_text )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
00696 (char*)attData );
00697 if( success )
00698 MB_SET_ERR( MB_FAILURE, "Failed to read byte/char data for attribute " << attIt->second.attName );
00699 ssAtt << "char;";
00700 break;
00701 case NC_SHORT:
00702 sz = attIt->second.attLen * sizeof( short );
00703 attData = (short*)malloc( sz );
00704 success = NCFUNC( get_att_short )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
00705 (short*)attData );
00706 if( success )
00707 MB_SET_ERR( MB_FAILURE, "Failed to read short data for attribute " << attIt->second.attName );
00708 ssAtt << "short;";
00709 break;
00710 case NC_INT:
00711 sz = attIt->second.attLen * sizeof( int );
00712 attData = (int*)malloc( sz );
00713 success = NCFUNC( get_att_int )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
00714 (int*)attData );
00715 if( success )
00716 MB_SET_ERR( MB_FAILURE, "Failed to read int data for attribute " << attIt->second.attName );
00717 ssAtt << "int;";
00718 break;
00719 case NC_FLOAT:
00720 sz = attIt->second.attLen * sizeof( float );
00721 attData = (float*)malloc( sz );
00722 success = NCFUNC( get_att_float )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
00723 (float*)attData );
00724 if( success )
00725 MB_SET_ERR( MB_FAILURE, "Failed to read float data for attribute " << attIt->second.attName );
00726 ssAtt << "float;";
00727 break;
00728 case NC_DOUBLE:
00729 sz = attIt->second.attLen * sizeof( double );
00730 attData = (double*)malloc( sz );
00731 success = NCFUNC( get_att_double )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
00732 (double*)attData );
00733 if( success )
00734 MB_SET_ERR( MB_FAILURE, "Failed to read double data for attribute " << attIt->second.attName );
00735 ssAtt << "double;";
00736 break;
00737 default:
00738 MB_SET_ERR( MB_FAILURE, "Unexpected data type for attribute " << attIt->second.attName );
00739 }
00740 char* tmpc = (char*)attData;
00741 for( unsigned int counter = 0; counter != sz; ++counter )
00742 ssAtt << tmpc[counter];
00743 free( attData );
00744 ssAtt << ';';
00745 attLen.push_back( ssAtt.str().size() - 1 );
00746 }
00747 attVal = ssAtt.str();
00748
00749 return MB_SUCCESS;
00750 }
00751
00752 ErrorCode NCHelper::create_dummy_variables()
00753 {
00754 Interface*& mbImpl = _readNC->mbImpl;
00755 std::vector< std::string >& dimNames = _readNC->dimNames;
00756 std::vector< int >& dimLens = _readNC->dimLens;
00757 std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
00758 DebugOutput& dbgOut = _readNC->dbgOut;
00759
00760 // Hack: look at all dimensions, and see if we have one that does not appear in the list of
00761 // varInfo names Right now, candidates are from unstructured meshes, such as ncol (HOMME) and
00762 // nCells (MPAS) For each of them, create a dummy coordinate variable with a sparse tag to store
00763 // the dimension length
00764 for( unsigned int i = 0; i < dimNames.size(); i++ )
00765 {
00766 // If there is a variable with this dimension name, skip
00767 if( varInfo.find( dimNames[i] ) != varInfo.end() ) continue;
00768
00769 // Create a dummy coordinate variable
00770 int sizeTotalVar = varInfo.size();
00771 std::string var_name( dimNames[i] );
00772 ReadNC::VarData& data = varInfo[var_name];
00773 data.varName = var_name;
00774 data.varId = sizeTotalVar;
00775 data.varTags.resize( 1, 0 );
00776 data.varDataType = NC_INT;
00777 data.varDims.resize( 1 );
00778 data.varDims[0] = (int)i;
00779 data.numAtts = 0;
00780 data.entLoc = ReadNC::ENTLOCSET;
00781 dummyVarNames.insert( var_name );
00782 dbgOut.tprintf( 2, "Dummy coordinate variable created for dimension %s\n", var_name.c_str() );
00783
00784 // Create a corresponding sparse tag
00785 Tag tagh;
00786 ErrorCode rval = mbImpl->tag_get_handle( var_name.c_str(), 0, MB_TYPE_INTEGER, tagh,
00787 MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag for dummy coordinate variable " << var_name );
00788
00789 // Tag value is the dimension length
00790 const void* ptr = &dimLens[i];
00791 // Tag size is 1
00792 int size = 1;
00793 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 );
00794
00795 dbgOut.tprintf( 2, "Sparse tag created for dimension %s\n", var_name.c_str() );
00796 }
00797
00798 return MB_SUCCESS;
00799 }
00800
00801 ErrorCode NCHelper::read_variables_to_set_allocate( std::vector< ReadNC::VarData >& vdatas,
00802 std::vector< int >& tstep_nums )
00803 {
00804 std::vector< int >& dimLens = _readNC->dimLens;
00805 DebugOutput& dbgOut = _readNC->dbgOut;
00806
00807 ErrorCode rval = MB_SUCCESS;
00808
00809 for( unsigned int i = 0; i < vdatas.size(); i++ )
00810 {
00811 // Set up readStarts and readCounts
00812 if( vdatas[i].has_tsteps )
00813 {
00814 // First: time
00815 vdatas[i].readStarts.push_back( 0 ); // This value is timestep dependent, will be set later
00816 vdatas[i].readCounts.push_back( 1 );
00817
00818 // Next: other dimensions
00819 for( unsigned int idx = 1; idx != vdatas[i].varDims.size(); idx++ )
00820 {
00821 vdatas[i].readStarts.push_back( 0 );
00822 vdatas[i].readCounts.push_back( dimLens[vdatas[i].varDims[idx]] );
00823 }
00824 }
00825 else
00826 {
00827 if( vdatas[i].varDims.empty() )
00828 {
00829 // Scalar variable
00830 vdatas[i].readStarts.push_back( 0 );
00831 vdatas[i].readCounts.push_back( 1 );
00832 }
00833 else
00834 {
00835 for( unsigned int idx = 0; idx != vdatas[i].varDims.size(); idx++ )
00836 {
00837 vdatas[i].readStarts.push_back( 0 );
00838 vdatas[i].readCounts.push_back( dimLens[vdatas[i].varDims[idx]] );
00839 }
00840 }
00841 }
00842
00843 // Get variable size
00844 vdatas[i].sz = 1;
00845 for( std::size_t idx = 0; idx != vdatas[i].readCounts.size(); idx++ )
00846 vdatas[i].sz *= vdatas[i].readCounts[idx];
00847
00848 // Note, for set variables without timesteps, loop one time and then break
00849 for( unsigned int t = 0; t < tstep_nums.size(); t++ )
00850 {
00851 dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
00852
00853 if( tstep_nums[t] >= dimLens[tDim] )
00854 {
00855 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] );
00856 }
00857
00858 // Get the tag to read into
00859 if( !vdatas[i].varTags[t] )
00860 {
00861 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 );
00862 }
00863
00864 switch( vdatas[i].varDataType )
00865 {
00866 case NC_BYTE:
00867 case NC_CHAR:
00868 vdatas[i].varDatas[t] = new char[vdatas[i].sz];
00869 break;
00870 case NC_SHORT:
00871 case NC_INT:
00872 vdatas[i].varDatas[t] = new int[vdatas[i].sz];
00873 break;
00874 case NC_FLOAT:
00875 case NC_DOUBLE:
00876 vdatas[i].varDatas[t] = new double[vdatas[i].sz];
00877 break;
00878 default:
00879 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
00880 }
00881
00882 // Loop continues only for set variables with timesteps, e.g. xtime(Time) or xtime(Time,
00883 // StrLen)
00884 if( !vdatas[i].has_tsteps ) break;
00885 }
00886 }
00887
00888 return rval;
00889 }
00890
00891 ErrorCode ScdNCHelper::check_existing_mesh()
00892 {
00893 Interface*& mbImpl = _readNC->mbImpl;
00894
00895 // Get the number of vertices
00896 int num_verts;
00897 ErrorCode rval = mbImpl->get_number_entities_by_dimension( _fileSet, 0, num_verts );MB_CHK_SET_ERR( rval, "Trouble getting number of vertices" );
00898
00899 /*
00900 // Check against parameters
00901 // When ghosting is used, this check might fail (to be updated later)
00902 if (num_verts > 0) {
00903 int expected_verts = (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1) * (-1 == lDims[2] ?
00904 1 : lDims[5] - lDims[2] + 1); if (num_verts != expected_verts) { MB_SET_ERR(MB_FAILURE, "Number
00905 of vertices doesn't match");
00906 }
00907 }
00908 */
00909
00910 // Check the number of elements too
00911 int num_elems;
00912 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" );
00913
00914 /*
00915 // Check against parameters
00916 // When ghosting is used, this check might fail (to be updated later)
00917 if (num_elems > 0) {
00918 int expected_elems = (lCDims[3] - lCDims[0] + 1) * (lCDims[4] - lCDims[1] + 1) * (-1 ==
00919 lCDims[2] ? 1 : (lCDims[5] - lCDims[2] + 1)); if (num_elems != expected_elems) {
00920 MB_SET_ERR(MB_FAILURE, "Number of elements doesn't match");
00921 }
00922 }
00923 */
00924
00925 return MB_SUCCESS;
00926 }
00927
00928 ErrorCode ScdNCHelper::create_mesh( Range& faces )
00929 {
00930 Interface*& mbImpl = _readNC->mbImpl;
00931 Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
00932 const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
00933 DebugOutput& dbgOut = _readNC->dbgOut;
00934 ScdInterface* scdi = _readNC->scdi;
00935 ScdParData& parData = _readNC->parData;
00936
00937 Range tmp_range;
00938 ScdBox* scd_box;
00939
00940 ErrorCode rval =
00941 scdi->construct_box( HomCoord( lDims[0], lDims[1], lDims[2], 1 ), HomCoord( lDims[3], lDims[4], lDims[5], 1 ),
00942 NULL, 0, scd_box, locallyPeriodic, &parData, true );MB_CHK_SET_ERR( rval, "Trouble creating scd vertex sequence" );
00943
00944 // Add verts to tmp_range first, so we can duplicate global ids in vertex ids
00945 tmp_range.insert( scd_box->start_vertex(), scd_box->start_vertex() + scd_box->num_vertices() - 1 );
00946
00947 if( mpFileIdTag )
00948 {
00949 int count;
00950 void* data;
00951 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" );
00952 assert( count == scd_box->num_vertices() );
00953 int* fid_data = (int*)data;
00954 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" );
00955 assert( count == scd_box->num_vertices() );
00956 int* gid_data = (int*)data;
00957 for( int i = 0; i < count; i++ )
00958 fid_data[i] = gid_data[i];
00959 }
00960
00961 // Then add box set and elements to the range, then to the file set
00962 tmp_range.insert( scd_box->start_element(), scd_box->start_element() + scd_box->num_elements() - 1 );
00963 tmp_range.insert( scd_box->box_set() );
00964 rval = mbImpl->add_entities( _fileSet, tmp_range );MB_CHK_SET_ERR( rval, "Couldn't add new vertices to current file set" );
00965
00966 dbgOut.tprintf( 1, "scdbox %d quads, %d vertices\n", scd_box->num_elements(), scd_box->num_vertices() );
00967
00968 // Set the vertex coordinates
00969 double *xc, *yc, *zc;
00970 rval = scd_box->get_coordinate_arrays( xc, yc, zc );MB_CHK_SET_ERR( rval, "Couldn't get vertex coordinate arrays" );
00971
00972 int i, j, k, il, jl, kl;
00973 int dil = lDims[3] - lDims[0] + 1;
00974 int djl = lDims[4] - lDims[1] + 1;
00975 assert( dil == (int)ilVals.size() && djl == (int)jlVals.size() &&
00976 ( -1 == lDims[2] || lDims[5] - lDims[2] + 1 == (int)levVals.size() ) );
00977
00978 for( kl = lDims[2]; kl <= lDims[5]; kl++ )
00979 {
00980 k = kl - lDims[2];
00981 for( jl = lDims[1]; jl <= lDims[4]; jl++ )
00982 {
00983 j = jl - lDims[1];
00984 for( il = lDims[0]; il <= lDims[3]; il++ )
00985 {
00986 i = il - lDims[0];
00987 unsigned int pos = i + j * dil + k * dil * djl;
00988 xc[pos] = ilVals[i];
00989 yc[pos] = jlVals[j];
00990 zc[pos] = ( -1 == lDims[2] ? 0.0 : levVals[k] );
00991 }
00992 }
00993 }
00994
00995 #ifndef NDEBUG
00996 int num_verts =
00997 ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) * ( -1 == lDims[2] ? 1 : lDims[5] - lDims[2] + 1 );
00998 std::vector< int > gids( num_verts );
00999 Range verts( scd_box->start_vertex(), scd_box->start_vertex() + scd_box->num_vertices() - 1 );
01000 rval = mbImpl->tag_get_data( mGlobalIdTag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting local gid values of vertices" );
01001 int vmin = *( std::min_element( gids.begin(), gids.end() ) ),
01002 vmax = *( std::max_element( gids.begin(), gids.end() ) );
01003 dbgOut.tprintf( 1, "Vertex gids %d-%d\n", vmin, vmax );
01004 #endif
01005
01006 // Add elements to the range passed in
01007 faces.insert( scd_box->start_element(), scd_box->start_element() + scd_box->num_elements() - 1 );
01008
01009 if( 2 <= dbgOut.get_verbosity() )
01010 {
01011 assert( scd_box->boundary_complete() );
01012 EntityHandle dum_ent = scd_box->start_element();
01013 rval = mbImpl->list_entities( &dum_ent, 1 );MB_CHK_SET_ERR( rval, "Trouble listing first hex" );
01014
01015 std::vector< EntityHandle > connect;
01016 rval = mbImpl->get_connectivity( &dum_ent, 1, connect );MB_CHK_SET_ERR( rval, "Trouble getting connectivity" );
01017
01018 rval = mbImpl->list_entities( &connect[0], connect.size() );MB_CHK_SET_ERR( rval, "Trouble listing element connectivity" );
01019 }
01020
01021 Range edges;
01022 mbImpl->get_adjacencies( faces, 1, true, edges, Interface::UNION );
01023
01024 // Create COORDS tag for quads
01025 rval = create_quad_coordinate_tag();MB_CHK_SET_ERR( rval, "Trouble creating COORDS tag for quads" );
01026
01027 return MB_SUCCESS;
01028 }
01029
01030 ErrorCode ScdNCHelper::read_variables( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
01031 {
01032 std::vector< ReadNC::VarData > vdatas;
01033 std::vector< ReadNC::VarData > vsetdatas;
01034
01035 ErrorCode rval = read_variables_setup( var_names, tstep_nums, vdatas, vsetdatas );MB_CHK_SET_ERR( rval, "Trouble setting up to read variables" );
01036
01037 if( !vsetdatas.empty() )
01038 {
01039 rval = read_variables_to_set( vsetdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to set" );
01040 }
01041
01042 if( !vdatas.empty() )
01043 {
01044 rval = read_scd_variables_to_nonset( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to verts/edges/faces" );
01045 }
01046
01047 return MB_SUCCESS;
01048 }
01049
01050 ErrorCode ScdNCHelper::read_scd_variables_to_nonset_allocate( std::vector< ReadNC::VarData >& vdatas,
01051 std::vector< int >& tstep_nums )
01052 {
01053 Interface*& mbImpl = _readNC->mbImpl;
01054 std::vector< int >& dimLens = _readNC->dimLens;
01055 DebugOutput& dbgOut = _readNC->dbgOut;
01056
01057 Range* range = NULL;
01058
01059 // Get vertices
01060 Range verts;
01061 ErrorCode rval = mbImpl->get_entities_by_dimension( _fileSet, 0, verts );MB_CHK_SET_ERR( rval, "Trouble getting vertices in current file set" );
01062 assert( "Should only have a single vertex subrange, since they were read in one shot" && verts.psize() == 1 );
01063
01064 Range edges;
01065 rval = mbImpl->get_entities_by_dimension( _fileSet, 1, edges );MB_CHK_SET_ERR( rval, "Trouble getting edges in current file set" );
01066
01067 // Get faces
01068 Range faces;
01069 rval = mbImpl->get_entities_by_dimension( _fileSet, 2, faces );MB_CHK_SET_ERR( rval, "Trouble getting faces in current file set" );
01070 assert( "Should only have a single face subrange, since they were read in one shot" && faces.psize() == 1 );
01071
01072 #ifdef MOAB_HAVE_MPI
01073 moab::Range faces_owned;
01074 bool& isParallel = _readNC->isParallel;
01075 if( isParallel )
01076 {
01077 ParallelComm*& myPcomm = _readNC->myPcomm;
01078 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" );
01079 }
01080 else
01081 faces_owned = faces; // Not running in parallel, but still with MPI
01082 #endif
01083
01084 for( unsigned int i = 0; i < vdatas.size(); i++ )
01085 {
01086 // Support non-set variables with 4 dimensions like (time, lev, lat, lon)
01087 assert( 4 == vdatas[i].varDims.size() );
01088
01089 // For a non-set variable, time should be the first dimension
01090 assert( tDim == vdatas[i].varDims[0] );
01091
01092 // Set up readStarts and readCounts
01093 vdatas[i].readStarts.resize( 4 );
01094 vdatas[i].readCounts.resize( 4 );
01095
01096 // First: time
01097 vdatas[i].readStarts[0] = 0; // This value is timestep dependent, will be set later
01098 vdatas[i].readCounts[0] = 1;
01099
01100 // Next: lev
01101 vdatas[i].readStarts[1] = 0;
01102 vdatas[i].readCounts[1] = vdatas[i].numLev;
01103
01104 // Finally: lat (or slat) and lon (or slon)
01105 switch( vdatas[i].entLoc )
01106 {
01107 case ReadNC::ENTLOCVERT:
01108 // Vertices
01109 vdatas[i].readStarts[2] = lDims[1];
01110 vdatas[i].readCounts[2] = lDims[4] - lDims[1] + 1;
01111 vdatas[i].readStarts[3] = lDims[0];
01112 vdatas[i].readCounts[3] = lDims[3] - lDims[0] + 1;
01113 range = &verts;
01114 break;
01115 case ReadNC::ENTLOCNSEDGE:
01116 case ReadNC::ENTLOCEWEDGE:
01117 case ReadNC::ENTLOCEDGE:
01118 // Not implemented yet, set a global error
01119 MB_SET_GLB_ERR( MB_NOT_IMPLEMENTED, "Reading edge data is not implemented yet" );
01120 case ReadNC::ENTLOCFACE:
01121 // Faces
01122 vdatas[i].readStarts[2] = lCDims[1];
01123 vdatas[i].readCounts[2] = lCDims[4] - lCDims[1] + 1;
01124 vdatas[i].readStarts[3] = lCDims[0];
01125 vdatas[i].readCounts[3] = lCDims[3] - lCDims[0] + 1;
01126 #ifdef MOAB_HAVE_MPI
01127 range = &faces_owned;
01128 #else
01129 range = &faces;
01130 #endif
01131 break;
01132 default:
01133 MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << vdatas[i].varName );
01134 }
01135
01136 for( unsigned int t = 0; t < tstep_nums.size(); t++ )
01137 {
01138 dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
01139
01140 if( tstep_nums[t] >= dimLens[tDim] )
01141 {
01142 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] );
01143 }
01144
01145 // Get the tag to read into
01146 if( !vdatas[i].varTags[t] )
01147 {
01148 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 );
01149 }
01150
01151 // Get ptr to tag space
01152 void* data;
01153 int count;
01154 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 );
01155 assert( (unsigned)count == range->size() );
01156 vdatas[i].varDatas[t] = data;
01157 }
01158
01159 // Get variable size
01160 vdatas[i].sz = 1;
01161 for( std::size_t idx = 0; idx != vdatas[i].readCounts.size(); idx++ )
01162 vdatas[i].sz *= vdatas[i].readCounts[idx];
01163 }
01164
01165 return rval;
01166 }
01167
01168 ErrorCode ScdNCHelper::read_scd_variables_to_nonset( std::vector< ReadNC::VarData >& vdatas,
01169 std::vector< int >& tstep_nums )
01170 {
01171 DebugOutput& dbgOut = _readNC->dbgOut;
01172
01173 ErrorCode rval = read_scd_variables_to_nonset_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read non-set variables" );
01174
01175 // Finally, read into that space
01176 int success;
01177 for( unsigned int i = 0; i < vdatas.size(); i++ )
01178 {
01179 std::size_t sz = vdatas[i].sz;
01180
01181 // A typical supported variable: float T(time, lev, lat, lon)
01182 // For tag values, need transpose (lev, lat, lon) to (lat, lon, lev)
01183 size_t ni = vdatas[i].readCounts[3]; // lon or slon
01184 size_t nj = vdatas[i].readCounts[2]; // lat or slat
01185 size_t nk = vdatas[i].readCounts[1]; // lev
01186
01187 for( unsigned int t = 0; t < tstep_nums.size(); t++ )
01188 {
01189 // Tag data for this timestep
01190 void* data = vdatas[i].varDatas[t];
01191
01192 // Set readStart for each timestep along time dimension
01193 vdatas[i].readStarts[0] = tstep_nums[t];
01194
01195 switch( vdatas[i].varDataType )
01196 {
01197 case NC_BYTE:
01198 case NC_CHAR: {
01199 std::vector< char > tmpchardata( sz );
01200 success = NCFUNCAG( _vara_text )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
01201 &vdatas[i].readCounts[0], &tmpchardata[0] );
01202 if( success )
01203 MB_SET_ERR( MB_FAILURE, "Failed to read byte/char data for variable " << vdatas[i].varName );
01204 if( vdatas[i].numLev > 1 )
01205 // Transpose (lev, lat, lon) to (lat, lon, lev)
01206 kji_to_jik( ni, nj, nk, data, &tmpchardata[0] );
01207 else
01208 {
01209 for( std::size_t idx = 0; idx != tmpchardata.size(); idx++ )
01210 ( (char*)data )[idx] = tmpchardata[idx];
01211 }
01212 break;
01213 }
01214 case NC_SHORT:
01215 case NC_INT: {
01216 std::vector< int > tmpintdata( sz );
01217 success = NCFUNCAG( _vara_int )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
01218 &vdatas[i].readCounts[0], &tmpintdata[0] );
01219 if( success )
01220 MB_SET_ERR( MB_FAILURE, "Failed to read short/int data for variable " << vdatas[i].varName );
01221 if( vdatas[i].numLev > 1 )
01222 // Transpose (lev, lat, lon) to (lat, lon, lev)
01223 kji_to_jik( ni, nj, nk, data, &tmpintdata[0] );
01224 else
01225 {
01226 for( std::size_t idx = 0; idx != tmpintdata.size(); idx++ )
01227 ( (int*)data )[idx] = tmpintdata[idx];
01228 }
01229 break;
01230 }
01231 case NC_FLOAT:
01232 case NC_DOUBLE: {
01233 std::vector< double > tmpdoubledata( sz );
01234 success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
01235 &vdatas[i].readCounts[0], &tmpdoubledata[0] );
01236 if( success )
01237 MB_SET_ERR( MB_FAILURE, "Failed to read float/double data for variable " << vdatas[i].varName );
01238 if( vdatas[i].numLev > 1 )
01239 // Transpose (lev, lat, lon) to (lat, lon, lev)
01240 kji_to_jik( ni, nj, nk, data, &tmpdoubledata[0] );
01241 else
01242 {
01243 for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
01244 ( (double*)data )[idx] = tmpdoubledata[idx];
01245 }
01246 break;
01247 }
01248 default:
01249 MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
01250 }
01251 }
01252 }
01253
01254 // Debug output, if requested
01255 if( 1 == dbgOut.get_verbosity() )
01256 {
01257 dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
01258 for( unsigned int i = 1; i < vdatas.size(); i++ )
01259 dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
01260 dbgOut.tprintf( 1, "\n" );
01261 }
01262
01263 return rval;
01264 }
01265
01266 ErrorCode ScdNCHelper::create_quad_coordinate_tag()
01267 {
01268 Interface*& mbImpl = _readNC->mbImpl;
01269
01270 Range ents;
01271 ErrorCode rval = mbImpl->get_entities_by_type( _fileSet, moab::MBQUAD, ents );MB_CHK_SET_ERR( rval, "Trouble getting quads" );
01272
01273 std::size_t numOwnedEnts = 0;
01274 #ifdef MOAB_HAVE_MPI
01275 Range ents_owned;
01276 bool& isParallel = _readNC->isParallel;
01277 if( isParallel )
01278 {
01279 ParallelComm*& myPcomm = _readNC->myPcomm;
01280 rval = myPcomm->filter_pstatus( ents, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &ents_owned );MB_CHK_SET_ERR( rval, "Trouble getting owned quads" );
01281 numOwnedEnts = ents_owned.size();
01282 }
01283 else
01284 {
01285 numOwnedEnts = ents.size();
01286 ents_owned = ents;
01287 }
01288 #else
01289 numOwnedEnts = ents.size();
01290 #endif
01291
01292 if( numOwnedEnts == 0 ) return MB_SUCCESS;
01293
01294 assert( numOwnedEnts == ilCVals.size() * jlCVals.size() );
01295 std::vector< double > coords( numOwnedEnts * 3 );
01296 std::size_t pos = 0;
01297 for( std::size_t j = 0; j != jlCVals.size(); ++j )
01298 {
01299 for( std::size_t i = 0; i != ilCVals.size(); ++i )
01300 {
01301 pos = j * ilCVals.size() * 3 + i * 3;
01302 coords[pos] = ilCVals[i];
01303 coords[pos + 1] = jlCVals[j];
01304 coords[pos + 2] = 0.0;
01305 }
01306 }
01307 std::string tag_name = "COORDS";
01308 Tag tagh = 0;
01309 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" );
01310
01311 void* data;
01312 int count;
01313 #ifdef MOAB_HAVE_MPI
01314 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" );
01315 #else
01316 rval = mbImpl->tag_iterate( tagh, ents.begin(), ents.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate COORDS tag on quads" );
01317 #endif
01318 assert( count == (int)numOwnedEnts );
01319 double* quad_data = (double*)data;
01320 std::copy( coords.begin(), coords.end(), quad_data );
01321
01322 return MB_SUCCESS;
01323 }
01324
01325 ErrorCode UcdNCHelper::read_variables( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
01326 {
01327 std::vector< ReadNC::VarData > vdatas;
01328 std::vector< ReadNC::VarData > vsetdatas;
01329
01330 ErrorCode rval = read_variables_setup( var_names, tstep_nums, vdatas, vsetdatas );MB_CHK_SET_ERR( rval, "Trouble setting up to read variables" );
01331
01332 if( !vsetdatas.empty() )
01333 {
01334 rval = read_variables_to_set( vsetdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to set" );
01335 }
01336
01337 if( !vdatas.empty() )
01338 {
01339 #ifdef MOAB_HAVE_PNETCDF
01340 // With pnetcdf support, we will use async read
01341 rval = read_ucd_variables_to_nonset_async( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to verts/edges/faces" );
01342 #else
01343 // Without pnetcdf support, we will use old read
01344 rval = read_ucd_variables_to_nonset( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to verts/edges/faces" );
01345 #endif
01346 }
01347
01348 return MB_SUCCESS;
01349 }
01350
01351 static double tolerance = 1.e-12;
01352 // Used to put points in an STL tree-based container
01353 bool NCHelper::Node3D::operator<( const NCHelper::Node3D& other ) const
01354 {
01355 if( coords[0] <= other.coords[0] - tolerance )
01356 return true;
01357 else if( coords[0] >= other.coords[0] + tolerance )
01358 return false;
01359 if( coords[1] <= other.coords[1] - tolerance )
01360 return true;
01361 else if( coords[1] >= other.coords[1] + tolerance )
01362 return false;
01363 if( coords[2] <= other.coords[2] - tolerance )
01364 return true;
01365 else if( coords[0] >= other.coords[0] + tolerance )
01366 return false;
01367
01368 return false;
01369 }
01370 /*bool NCHelper::Node3D::operator == ( const NCHelper::Node3D& other ) const
01371 return ( !(*this < other) && ! ( other < *this) ) ;*/
01372 } // namespace moab
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.