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