MOAB: Mesh Oriented datABase  (version 5.2.1)
ReadNC.cpp
Go to the documentation of this file.
00001 #include "ReadNC.hpp"
00002 #include "NCHelper.hpp"
00003 
00004 #include "moab/ReadUtilIface.hpp"
00005 #include "MBTagConventions.hpp"
00006 #include "moab/FileOptions.hpp"
00007 
00008 namespace moab
00009 {
00010 
00011 ReaderIface* ReadNC::factory( Interface* iface )
00012 {
00013     return new ReadNC( iface );
00014 }
00015 
00016 ReadNC::ReadNC( Interface* impl )
00017     : mbImpl( impl ), fileId( -1 ), mGlobalIdTag( 0 ), mpFileIdTag( NULL ), dbgOut( stderr ), isParallel( false ),
00018       partMethod( ScdParData::ALLJORKORI ), scdi( NULL ),
00019 #ifdef MOAB_HAVE_MPI
00020       myPcomm( NULL ),
00021 #endif
00022       noMesh( false ), noVars( false ), spectralMesh( false ), noMixedElements( false ), noEdges( false ),
00023       gatherSetRank( -1 ), tStepBase( -1 ), trivialPartitionShift( 0 ), myHelper( NULL )
00024 {
00025     assert( impl != NULL );
00026     impl->query_interface( readMeshIface );
00027 }
00028 
00029 ReadNC::~ReadNC()
00030 {
00031     mbImpl->release_interface( readMeshIface );
00032     if( myHelper != NULL ) delete myHelper;
00033 }
00034 
00035 ErrorCode ReadNC::load_file( const char* file_name, const EntityHandle* file_set, const FileOptions& opts,
00036                              const ReaderIface::SubsetList* /*subset_list*/, const Tag* file_id_tag )
00037 {
00038     // See if opts has variable(s) specified
00039     std::vector< std::string > var_names;
00040     std::vector< int > tstep_nums;
00041     std::vector< double > tstep_vals;
00042 
00043     // Get and cache predefined tag handles
00044     mGlobalIdTag = mbImpl->globalId_tag();
00045     // Store the pointer to the tag; if not null, set when global id tag
00046     // is set too, with the same data, duplicated
00047     mpFileIdTag = file_id_tag;
00048 
00049     ErrorCode rval = parse_options( opts, var_names, tstep_nums, tstep_vals );MB_CHK_SET_ERR( rval, "Trouble parsing option string" );
00050 
00051     // Open the file
00052     dbgOut.tprintf( 1, "Opening file %s\n", file_name );
00053     fileName = std::string( file_name );
00054     int success;
00055 
00056 #ifdef MOAB_HAVE_PNETCDF
00057     if( isParallel )
00058         success = NCFUNC( open )( myPcomm->proc_config().proc_comm(), file_name, 0, MPI_INFO_NULL, &fileId );
00059     else
00060         success = NCFUNC( open )( MPI_COMM_SELF, file_name, 0, MPI_INFO_NULL, &fileId );
00061 #else
00062     success = NCFUNC( open )( file_name, 0, &fileId );
00063 #endif
00064     if( success ) MB_SET_ERR( MB_FAILURE, "Trouble opening file " << file_name );
00065 
00066     // Read the header (num dimensions, dimensions, num variables, global attribs)
00067     rval = read_header();MB_CHK_SET_ERR( rval, "Trouble reading file header" );
00068 
00069     // Make sure there's a file set to put things in
00070     EntityHandle tmp_set;
00071     if( noMesh && !file_set ) { MB_SET_ERR( MB_FAILURE, "NOMESH option requires non-NULL file set on input" ); }
00072     else if( !file_set || ( file_set && *file_set == 0 ) )
00073     {
00074         rval = mbImpl->create_meshset( MESHSET_SET, tmp_set );MB_CHK_SET_ERR( rval, "Trouble creating file set" );
00075     }
00076     else
00077         tmp_set = *file_set;
00078 
00079     // Get the scd interface
00080     scdi = NULL;
00081     rval = mbImpl->query_interface( scdi );
00082     if( NULL == scdi ) return MB_FAILURE;
00083 
00084     if( NULL != myHelper ) delete myHelper;
00085 
00086     // Get appropriate NC helper instance based on information read from the header
00087     myHelper = NCHelper::get_nc_helper( this, fileId, opts, tmp_set );
00088     if( NULL == myHelper ) { MB_SET_ERR( MB_FAILURE, "Failed to get NCHelper class instance" ); }
00089 
00090     // Initialize mesh values
00091     rval = myHelper->init_mesh_vals();MB_CHK_SET_ERR( rval, "Trouble initializing mesh values" );
00092 
00093     // Check existing mesh from last read
00094     if( noMesh && !noVars )
00095     {
00096         rval = myHelper->check_existing_mesh();MB_CHK_SET_ERR( rval, "Trouble checking mesh from last read" );
00097     }
00098 
00099     // Create some conventional tags, e.g. __NUM_DIMS
00100     // For multiple reads to a specified file set, we assume a single file, or a series of
00101     // files with separated timesteps. Keep a flag on the file set to prevent conventional
00102     // tags from being created again on a second read
00103     Tag convTagsCreated = 0;
00104     int def_val         = 0;
00105     rval                = mbImpl->tag_get_handle( "__CONV_TAGS_CREATED", 1, MB_TYPE_INTEGER, convTagsCreated,
00106                                    MB_TAG_SPARSE | MB_TAG_CREAT, &def_val );MB_CHK_SET_ERR( rval, "Trouble getting _CONV_TAGS_CREATED tag" );
00107     int create_conv_tags_flag = 0;
00108     rval                      = mbImpl->tag_get_data( convTagsCreated, &tmp_set, 1, &create_conv_tags_flag );
00109     // The first read to the file set
00110     if( 0 == create_conv_tags_flag )
00111     {
00112         // Read dimensions (coordinate variables) by default to create tags like __<var_name>_DIMS
00113         // This is done only once (assume that all files read to the file set have the same
00114         // dimensions)
00115         rval = myHelper->read_variables( dimNames, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading dimensions" );
00116 
00117         rval = myHelper->create_conventional_tags( tstep_nums );MB_CHK_SET_ERR( rval, "Trouble creating NC conventional tags" );
00118 
00119         create_conv_tags_flag = 1;
00120         rval                  = mbImpl->tag_set_data( convTagsCreated, &tmp_set, 1, &create_conv_tags_flag );MB_CHK_SET_ERR( rval, "Trouble setting data to _CONV_TAGS_CREATED tag" );
00121     }
00122     // Another read to the file set
00123     else
00124     {
00125         if( tStepBase > -1 )
00126         {
00127             // If timesteps spread across files, merge time values read
00128             // from current file to existing time tag
00129             rval = myHelper->update_time_tag_vals();MB_CHK_SET_ERR( rval, "Trouble updating time tag values" );
00130         }
00131     }
00132 
00133     // Create mesh vertex/edge/face sequences
00134     Range faces;
00135     if( !noMesh )
00136     {
00137         rval = myHelper->create_mesh( faces );MB_CHK_SET_ERR( rval, "Trouble creating mesh" );
00138     }
00139 
00140     // Read specified variables onto grid
00141     if( !noVars )
00142     {
00143         if( var_names.empty() )
00144         {
00145             // If VARIABLE option is missing, read all variables
00146             rval = myHelper->read_variables( var_names, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading all variables" );
00147         }
00148         else
00149         {
00150             // Exclude dimensions that are read to the file set by default
00151             std::vector< std::string > non_dim_var_names;
00152             for( unsigned int i = 0; i < var_names.size(); i++ )
00153             {
00154                 if( std::find( dimNames.begin(), dimNames.end(), var_names[i] ) == dimNames.end() )
00155                     non_dim_var_names.push_back( var_names[i] );
00156             }
00157 
00158             if( !non_dim_var_names.empty() )
00159             {
00160                 rval = myHelper->read_variables( non_dim_var_names, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading specified variables" );
00161             }
00162         }
00163     }
00164 
00165 #ifdef MOAB_HAVE_MPI
00166     // Create partition set, and populate with elements
00167     if( isParallel )
00168     {
00169         EntityHandle partn_set;
00170         rval = mbImpl->create_meshset( MESHSET_SET, partn_set );MB_CHK_SET_ERR( rval, "Trouble creating partition set" );
00171 
00172         rval = mbImpl->add_entities( partn_set, faces );MB_CHK_SET_ERR( rval, "Couldn't add new faces to partition set" );
00173 
00174         Range verts;
00175         rval = mbImpl->get_connectivity( faces, verts );MB_CHK_SET_ERR( rval, "Couldn't get verts of faces" );
00176 
00177         rval = mbImpl->add_entities( partn_set, verts );MB_CHK_SET_ERR( rval, "Couldn't add new verts to partition set" );
00178 
00179         myPcomm->partition_sets().insert( partn_set );
00180 
00181         // Write partition tag name on partition set
00182         Tag part_tag = myPcomm->partition_tag();
00183         int dum_rank = myPcomm->proc_config().proc_rank();
00184         rval         = mbImpl->tag_set_data( part_tag, &partn_set, 1, &dum_rank );MB_CHK_SET_ERR( rval, "Trouble writing partition tag name on partition set" );
00185     }
00186 #endif
00187 
00188     mbImpl->release_interface( scdi );
00189     scdi = NULL;
00190 
00191     // Close the file
00192     success = NCFUNC( close )( fileId );
00193     if( success ) MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
00194 
00195     return MB_SUCCESS;
00196 }
00197 
00198 ErrorCode ReadNC::parse_options( const FileOptions& opts, std::vector< std::string >& var_names,
00199                                  std::vector< int >& tstep_nums, std::vector< double >& tstep_vals )
00200 {
00201     int tmpval;
00202     if( MB_SUCCESS == opts.get_int_option( "DEBUG_IO", 1, tmpval ) )
00203     {
00204         dbgOut.set_verbosity( tmpval );
00205         dbgOut.set_prefix( "NC " );
00206     }
00207 
00208     ErrorCode rval = opts.get_strs_option( "VARIABLE", var_names );
00209     if( MB_TYPE_OUT_OF_RANGE == rval )
00210         noVars = true;
00211     else
00212         noVars = false;
00213 
00214     opts.get_ints_option( "TIMESTEP", tstep_nums );
00215     opts.get_reals_option( "TIMEVAL", tstep_vals );
00216 
00217     rval = opts.get_null_option( "NOMESH" );
00218     if( MB_SUCCESS == rval ) noMesh = true;
00219 
00220     rval = opts.get_null_option( "SPECTRAL_MESH" );
00221     if( MB_SUCCESS == rval ) spectralMesh = true;
00222 
00223     rval = opts.get_null_option( "NO_MIXED_ELEMENTS" );
00224     if( MB_SUCCESS == rval ) noMixedElements = true;
00225 
00226     rval = opts.get_null_option( "NO_EDGES" );
00227     if( MB_SUCCESS == rval ) noEdges = true;
00228 
00229     if( 2 <= dbgOut.get_verbosity() )
00230     {
00231         if( !var_names.empty() )
00232         {
00233             std::cerr << "Variables requested: ";
00234             for( unsigned int i = 0; i < var_names.size(); i++ )
00235                 std::cerr << var_names[i];
00236             std::cerr << std::endl;
00237         }
00238 
00239         if( !tstep_nums.empty() )
00240         {
00241             std::cerr << "Timesteps requested: ";
00242             for( unsigned int i = 0; i < tstep_nums.size(); i++ )
00243                 std::cerr << tstep_nums[i];
00244             std::cerr << std::endl;
00245         }
00246 
00247         if( !tstep_vals.empty() )
00248         {
00249             std::cerr << "Time vals requested: ";
00250             for( unsigned int i = 0; i < tstep_vals.size(); i++ )
00251                 std::cerr << tstep_vals[i];
00252             std::cerr << std::endl;
00253         }
00254     }
00255 
00256     rval = opts.get_int_option( "GATHER_SET", 0, gatherSetRank );
00257     if( MB_TYPE_OUT_OF_RANGE == rval ) { MB_SET_ERR( rval, "Invalid value for GATHER_SET option" ); }
00258 
00259     rval = opts.get_int_option( "TIMESTEPBASE", 0, tStepBase );
00260     if( MB_TYPE_OUT_OF_RANGE == rval ) { MB_SET_ERR( rval, "Invalid value for TIMESTEPBASE option" ); }
00261 
00262     rval = opts.get_int_option( "TRIVIAL_PARTITION_SHIFT", 1, trivialPartitionShift );
00263     if( MB_TYPE_OUT_OF_RANGE == rval ) { MB_SET_ERR( rval, "Invalid value for TRIVIAL_PARTITION_SHIFT option" ); }
00264 
00265 #ifdef MOAB_HAVE_MPI
00266     isParallel = ( opts.match_option( "PARALLEL", "READ_PART" ) != MB_ENTITY_NOT_FOUND );
00267 
00268     if( !isParallel )
00269         // Return success here, since rval still has _NOT_FOUND from not finding option
00270         // in this case, myPcomm will be NULL, so it can never be used; always check for isParallel
00271         // before any use for myPcomm
00272         return MB_SUCCESS;
00273 
00274     int pcomm_no = 0;
00275     rval         = opts.get_int_option( "PARALLEL_COMM", pcomm_no );
00276     if( MB_TYPE_OUT_OF_RANGE == rval ) { MB_SET_ERR( rval, "Invalid value for PARALLEL_COMM option" ); }
00277     myPcomm = ParallelComm::get_pcomm( mbImpl, pcomm_no );
00278     if( 0 == myPcomm ) { myPcomm = new ParallelComm( mbImpl, MPI_COMM_WORLD ); }
00279     const int rank = myPcomm->proc_config().proc_rank();
00280     dbgOut.set_rank( rank );
00281 
00282     int dum;
00283     rval = opts.match_option( "PARTITION_METHOD", ScdParData::PartitionMethodNames, dum );
00284     if( MB_FAILURE == rval ) { MB_SET_ERR( rval, "Unknown partition method specified" ); }
00285     else if( MB_ENTITY_NOT_FOUND == rval )
00286         partMethod = ScdParData::ALLJORKORI;
00287     else
00288         partMethod = dum;
00289 #endif
00290 
00291     return MB_SUCCESS;
00292 }
00293 
00294 ErrorCode ReadNC::read_header()
00295 {
00296     dbgOut.tprint( 1, "Reading header...\n" );
00297 
00298     // Get the global attributes
00299     int numgatts;
00300     int success;
00301     success = NCFUNC( inq_natts )( fileId, &numgatts );
00302     if( success ) MB_SET_ERR( MB_FAILURE, "Couldn't get number of global attributes" );
00303 
00304     // Read attributes into globalAtts
00305     ErrorCode result = get_attributes( NC_GLOBAL, numgatts, globalAtts );MB_CHK_SET_ERR( result, "Trouble getting global attributes" );
00306     dbgOut.tprintf( 1, "Read %u attributes\n", (unsigned int)globalAtts.size() );
00307 
00308     // Read in dimensions into dimNames and dimLens
00309     result = get_dimensions( fileId, dimNames, dimLens );MB_CHK_SET_ERR( result, "Trouble getting dimensions" );
00310     dbgOut.tprintf( 1, "Read %u dimensions\n", (unsigned int)dimNames.size() );
00311 
00312     // Read in variables into varInfo
00313     result = get_variables();MB_CHK_SET_ERR( result, "Trouble getting variables" );
00314     dbgOut.tprintf( 1, "Read %u variables\n", (unsigned int)varInfo.size() );
00315 
00316     return MB_SUCCESS;
00317 }
00318 
00319 ErrorCode ReadNC::get_attributes( int var_id, int num_atts, std::map< std::string, AttData >& atts, const char* prefix )
00320 {
00321     char dum_name[120];
00322 
00323     for( int i = 0; i < num_atts; i++ )
00324     {
00325         // Get the name
00326         int success = NCFUNC( inq_attname )( fileId, var_id, i, dum_name );
00327         if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting attribute name" );
00328 
00329         AttData& data = atts[std::string( dum_name )];
00330         data.attName  = std::string( dum_name );
00331         success       = NCFUNC( inq_att )( fileId, var_id, dum_name, &data.attDataType, &data.attLen );
00332         if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting info for attribute " << data.attName );
00333         data.attVarId = var_id;
00334 
00335         dbgOut.tprintf( 2, "%sAttribute %s: length=%u, varId=%d, type=%d\n", ( prefix ? prefix : "" ),
00336                         data.attName.c_str(), (unsigned int)data.attLen, data.attVarId, data.attDataType );
00337     }
00338 
00339     return MB_SUCCESS;
00340 }
00341 
00342 ErrorCode ReadNC::get_dimensions( int file_id, std::vector< std::string >& dim_names, std::vector< int >& dim_lens )
00343 {
00344     // Get the number of dimensions
00345     int num_dims;
00346     int success = NCFUNC( inq_ndims )( file_id, &num_dims );
00347     if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of dimensions" );
00348 
00349     if( num_dims > NC_MAX_DIMS )
00350     {
00351         MB_SET_ERR( MB_FAILURE,
00352                     "ReadNC: File contains " << num_dims << " dims but NetCDF library supports only " << NC_MAX_DIMS );
00353     }
00354 
00355     char dim_name[NC_MAX_NAME + 1];
00356     NCDF_SIZE dim_len;
00357     dim_names.resize( num_dims );
00358     dim_lens.resize( num_dims );
00359 
00360     for( int i = 0; i < num_dims; i++ )
00361     {
00362         success = NCFUNC( inq_dim )( file_id, i, dim_name, &dim_len );
00363         if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting dimension info" );
00364 
00365         dim_names[i] = std::string( dim_name );
00366         dim_lens[i]  = dim_len;
00367 
00368         dbgOut.tprintf( 2, "Dimension %s, length=%u\n", dim_name, (unsigned int)dim_len );
00369     }
00370 
00371     return MB_SUCCESS;
00372 }
00373 
00374 ErrorCode ReadNC::get_variables()
00375 {
00376     // First cache the number of time steps
00377     std::vector< std::string >::iterator vit = std::find( dimNames.begin(), dimNames.end(), "time" );
00378     if( vit == dimNames.end() ) vit = std::find( dimNames.begin(), dimNames.end(), "t" );
00379 
00380     int ntimes = 0;
00381     if( vit != dimNames.end() ) ntimes = dimLens[vit - dimNames.begin()];
00382     if( !ntimes ) ntimes = 1;
00383 
00384     // Get the number of variables
00385     int num_vars;
00386     int success = NCFUNC( inq_nvars )( fileId, &num_vars );
00387     if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of variables" );
00388 
00389     if( num_vars > NC_MAX_VARS )
00390     {
00391         MB_SET_ERR( MB_FAILURE,
00392                     "ReadNC: File contains " << num_vars << " vars but NetCDF library supports only " << NC_MAX_VARS );
00393     }
00394 
00395     char var_name[NC_MAX_NAME + 1];
00396     int var_ndims;
00397 
00398     for( int i = 0; i < num_vars; i++ )
00399     {
00400         // Get the name first, so we can allocate a map iterate for this var
00401         success = NCFUNC( inq_varname )( fileId, i, var_name );
00402         if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting variable name" );
00403         VarData& data = varInfo[std::string( var_name )];
00404         data.varName  = std::string( var_name );
00405         data.varId    = i;
00406         data.varTags.resize( ntimes, 0 );
00407 
00408         // Get the data type
00409         success = NCFUNC( inq_vartype )( fileId, i, &data.varDataType );
00410         if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting data type for variable " << data.varName );
00411 
00412         // Get the number of dimensions, then the dimensions
00413         success = NCFUNC( inq_varndims )( fileId, i, &var_ndims );
00414         if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of dims for variable " << data.varName );
00415         data.varDims.resize( var_ndims );
00416 
00417         success = NCFUNC( inq_vardimid )( fileId, i, &data.varDims[0] );
00418         if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting dimensions for variable " << data.varName );
00419 
00420         // Finally, get the number of attributes, then the attributes
00421         success = NCFUNC( inq_varnatts )( fileId, i, &data.numAtts );
00422         if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of dims for variable " << data.varName );
00423 
00424         // Print debug info here so attribute info comes afterwards
00425         dbgOut.tprintf( 2, "Variable %s: Id=%d, numAtts=%d, datatype=%d, num_dims=%u\n", data.varName.c_str(),
00426                         data.varId, data.numAtts, data.varDataType, (unsigned int)data.varDims.size() );
00427 
00428         ErrorCode rval = get_attributes( i, data.numAtts, data.varAtts, "   " );MB_CHK_SET_ERR( rval, "Trouble getting attributes for variable " << data.varName );
00429     }
00430 
00431     return MB_SUCCESS;
00432 }
00433 
00434 ErrorCode ReadNC::read_tag_values( const char*, const char*, const FileOptions&, std::vector< int >&,
00435                                    const SubsetList* )
00436 {
00437     return MB_FAILURE;
00438 }
00439 
00440 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines