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