MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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 // Write partition tag name on partition set 00170 Tag part_tag = myPcomm->partition_tag(); 00171 int dum_rank = myPcomm->proc_config().proc_rank(); 00172 // the tmp_set is the file_set 00173 rval = mbImpl->tag_set_data( part_tag, &tmp_set, 1, &dum_rank );MB_CHK_SET_ERR( rval, "Trouble writing partition tag name on partition set" ); 00174 } 00175 #endif 00176 00177 mbImpl->release_interface( scdi ); 00178 scdi = NULL; 00179 00180 // Close the file 00181 success = NCFUNC( close )( fileId ); 00182 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble closing file" ); 00183 00184 return MB_SUCCESS; 00185 } 00186 00187 ErrorCode ReadNC::parse_options( const FileOptions& opts, std::vector< std::string >& var_names, 00188 std::vector< int >& tstep_nums, std::vector< double >& tstep_vals ) 00189 { 00190 int tmpval; 00191 if( MB_SUCCESS == opts.get_int_option( "DEBUG_IO", 1, tmpval ) ) 00192 { 00193 dbgOut.set_verbosity( tmpval ); 00194 dbgOut.set_prefix( "NC " ); 00195 } 00196 00197 ErrorCode rval = opts.get_strs_option( "VARIABLE", var_names ); 00198 if( MB_TYPE_OUT_OF_RANGE == rval ) 00199 noVars = true; 00200 else 00201 noVars = false; 00202 00203 opts.get_ints_option( "TIMESTEP", tstep_nums ); 00204 opts.get_reals_option( "TIMEVAL", tstep_vals ); 00205 00206 rval = opts.get_null_option( "NOMESH" ); 00207 if( MB_SUCCESS == rval ) noMesh = true; 00208 00209 rval = opts.get_null_option( "SPECTRAL_MESH" ); 00210 if( MB_SUCCESS == rval ) spectralMesh = true; 00211 00212 rval = opts.get_null_option( "NO_MIXED_ELEMENTS" ); 00213 if( MB_SUCCESS == rval ) noMixedElements = true; 00214 00215 rval = opts.get_null_option( "NO_EDGES" ); 00216 if( MB_SUCCESS == rval ) noEdges = true; 00217 00218 if( 2 <= dbgOut.get_verbosity() ) 00219 { 00220 if( !var_names.empty() ) 00221 { 00222 std::cerr << "Variables requested: "; 00223 for( unsigned int i = 0; i < var_names.size(); i++ ) 00224 std::cerr << var_names[i]; 00225 std::cerr << std::endl; 00226 } 00227 00228 if( !tstep_nums.empty() ) 00229 { 00230 std::cerr << "Timesteps requested: "; 00231 for( unsigned int i = 0; i < tstep_nums.size(); i++ ) 00232 std::cerr << tstep_nums[i]; 00233 std::cerr << std::endl; 00234 } 00235 00236 if( !tstep_vals.empty() ) 00237 { 00238 std::cerr << "Time vals requested: "; 00239 for( unsigned int i = 0; i < tstep_vals.size(); i++ ) 00240 std::cerr << tstep_vals[i]; 00241 std::cerr << std::endl; 00242 } 00243 } 00244 00245 rval = opts.get_int_option( "GATHER_SET", 0, gatherSetRank ); 00246 if( MB_TYPE_OUT_OF_RANGE == rval ) { MB_SET_ERR( rval, "Invalid value for GATHER_SET option" ); } 00247 00248 rval = opts.get_int_option( "TIMESTEPBASE", 0, tStepBase ); 00249 if( MB_TYPE_OUT_OF_RANGE == rval ) { MB_SET_ERR( rval, "Invalid value for TIMESTEPBASE option" ); } 00250 00251 rval = opts.get_int_option( "TRIVIAL_PARTITION_SHIFT", 1, trivialPartitionShift ); 00252 if( MB_TYPE_OUT_OF_RANGE == rval ) { MB_SET_ERR( rval, "Invalid value for TRIVIAL_PARTITION_SHIFT option" ); } 00253 00254 #ifdef MOAB_HAVE_MPI 00255 isParallel = ( opts.match_option( "PARALLEL", "READ_PART" ) != MB_ENTITY_NOT_FOUND ); 00256 00257 if( !isParallel ) 00258 // Return success here, since rval still has _NOT_FOUND from not finding option 00259 // in this case, myPcomm will be NULL, so it can never be used; always check for isParallel 00260 // before any use for myPcomm 00261 return MB_SUCCESS; 00262 00263 int pcomm_no = 0; 00264 rval = opts.get_int_option( "PARALLEL_COMM", pcomm_no ); 00265 if( MB_TYPE_OUT_OF_RANGE == rval ) { MB_SET_ERR( rval, "Invalid value for PARALLEL_COMM option" ); } 00266 myPcomm = ParallelComm::get_pcomm( mbImpl, pcomm_no ); 00267 if( 0 == myPcomm ) { myPcomm = new ParallelComm( mbImpl, MPI_COMM_WORLD ); } 00268 const int rank = myPcomm->proc_config().proc_rank(); 00269 dbgOut.set_rank( rank ); 00270 00271 int dum; 00272 rval = opts.match_option( "PARTITION_METHOD", ScdParData::PartitionMethodNames, dum ); 00273 if( MB_FAILURE == rval ) { MB_SET_ERR( rval, "Unknown partition method specified" ); } 00274 else if( MB_ENTITY_NOT_FOUND == rval ) 00275 partMethod = ScdParData::ALLJORKORI; 00276 else 00277 partMethod = dum; 00278 #endif 00279 00280 return MB_SUCCESS; 00281 } 00282 00283 ErrorCode ReadNC::read_header() 00284 { 00285 dbgOut.tprint( 1, "Reading header...\n" ); 00286 00287 // Get the global attributes 00288 int numgatts; 00289 int success; 00290 success = NCFUNC( inq_natts )( fileId, &numgatts ); 00291 if( success ) MB_SET_ERR( MB_FAILURE, "Couldn't get number of global attributes" ); 00292 00293 // Read attributes into globalAtts 00294 ErrorCode result = get_attributes( NC_GLOBAL, numgatts, globalAtts );MB_CHK_SET_ERR( result, "Trouble getting global attributes" ); 00295 dbgOut.tprintf( 1, "Read %u attributes\n", (unsigned int)globalAtts.size() ); 00296 00297 // Read in dimensions into dimNames and dimLens 00298 result = get_dimensions( fileId, dimNames, dimLens );MB_CHK_SET_ERR( result, "Trouble getting dimensions" ); 00299 dbgOut.tprintf( 1, "Read %u dimensions\n", (unsigned int)dimNames.size() ); 00300 00301 // Read in variables into varInfo 00302 result = get_variables();MB_CHK_SET_ERR( result, "Trouble getting variables" ); 00303 dbgOut.tprintf( 1, "Read %u variables\n", (unsigned int)varInfo.size() ); 00304 00305 return MB_SUCCESS; 00306 } 00307 00308 ErrorCode ReadNC::get_attributes( int var_id, int num_atts, std::map< std::string, AttData >& atts, const char* prefix ) 00309 { 00310 char dum_name[120]; 00311 00312 for( int i = 0; i < num_atts; i++ ) 00313 { 00314 // Get the name 00315 int success = NCFUNC( inq_attname )( fileId, var_id, i, dum_name ); 00316 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting attribute name" ); 00317 00318 AttData& data = atts[std::string( dum_name )]; 00319 data.attName = std::string( dum_name ); 00320 success = NCFUNC( inq_att )( fileId, var_id, dum_name, &data.attDataType, &data.attLen ); 00321 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting info for attribute " << data.attName ); 00322 data.attVarId = var_id; 00323 00324 dbgOut.tprintf( 2, "%sAttribute %s: length=%u, varId=%d, type=%d\n", ( prefix ? prefix : "" ), 00325 data.attName.c_str(), (unsigned int)data.attLen, data.attVarId, data.attDataType ); 00326 } 00327 00328 return MB_SUCCESS; 00329 } 00330 00331 ErrorCode ReadNC::get_dimensions( int file_id, std::vector< std::string >& dim_names, std::vector< int >& dim_lens ) 00332 { 00333 // Get the number of dimensions 00334 int num_dims; 00335 int success = NCFUNC( inq_ndims )( file_id, &num_dims ); 00336 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of dimensions" ); 00337 00338 if( num_dims > NC_MAX_DIMS ) 00339 { 00340 MB_SET_ERR( MB_FAILURE, 00341 "ReadNC: File contains " << num_dims << " dims but NetCDF library supports only " << NC_MAX_DIMS ); 00342 } 00343 00344 char dim_name[NC_MAX_NAME + 1]; 00345 NCDF_SIZE dim_len; 00346 dim_names.resize( num_dims ); 00347 dim_lens.resize( num_dims ); 00348 00349 for( int i = 0; i < num_dims; i++ ) 00350 { 00351 success = NCFUNC( inq_dim )( file_id, i, dim_name, &dim_len ); 00352 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting dimension info" ); 00353 00354 dim_names[i] = std::string( dim_name ); 00355 dim_lens[i] = dim_len; 00356 00357 dbgOut.tprintf( 2, "Dimension %s, length=%u\n", dim_name, (unsigned int)dim_len ); 00358 } 00359 00360 return MB_SUCCESS; 00361 } 00362 00363 ErrorCode ReadNC::get_variables() 00364 { 00365 // First cache the number of time steps 00366 std::vector< std::string >::iterator vit = std::find( dimNames.begin(), dimNames.end(), "time" ); 00367 if( vit == dimNames.end() ) vit = std::find( dimNames.begin(), dimNames.end(), "t" ); 00368 00369 int ntimes = 0; 00370 if( vit != dimNames.end() ) ntimes = dimLens[vit - dimNames.begin()]; 00371 if( !ntimes ) ntimes = 1; 00372 00373 // Get the number of variables 00374 int num_vars; 00375 int success = NCFUNC( inq_nvars )( fileId, &num_vars ); 00376 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of variables" ); 00377 00378 if( num_vars > NC_MAX_VARS ) 00379 { 00380 MB_SET_ERR( MB_FAILURE, 00381 "ReadNC: File contains " << num_vars << " vars but NetCDF library supports only " << NC_MAX_VARS ); 00382 } 00383 00384 char var_name[NC_MAX_NAME + 1]; 00385 int var_ndims; 00386 00387 for( int i = 0; i < num_vars; i++ ) 00388 { 00389 // Get the name first, so we can allocate a map iterate for this var 00390 success = NCFUNC( inq_varname )( fileId, i, var_name ); 00391 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting variable name" ); 00392 VarData& data = varInfo[std::string( var_name )]; 00393 data.varName = std::string( var_name ); 00394 data.varId = i; 00395 data.varTags.resize( ntimes, 0 ); 00396 00397 // Get the data type 00398 success = NCFUNC( inq_vartype )( fileId, i, &data.varDataType ); 00399 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting data type for variable " << data.varName ); 00400 00401 // Get the number of dimensions, then the dimensions 00402 success = NCFUNC( inq_varndims )( fileId, i, &var_ndims ); 00403 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of dims for variable " << data.varName ); 00404 data.varDims.resize( var_ndims ); 00405 00406 success = NCFUNC( inq_vardimid )( fileId, i, &data.varDims[0] ); 00407 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting dimensions for variable " << data.varName ); 00408 00409 // Finally, get the number of attributes, then the attributes 00410 success = NCFUNC( inq_varnatts )( fileId, i, &data.numAtts ); 00411 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble getting number of dims for variable " << data.varName ); 00412 00413 // Print debug info here so attribute info comes afterwards 00414 dbgOut.tprintf( 2, "Variable %s: Id=%d, numAtts=%d, datatype=%d, num_dims=%u\n", data.varName.c_str(), 00415 data.varId, data.numAtts, data.varDataType, (unsigned int)data.varDims.size() ); 00416 00417 ErrorCode rval = get_attributes( i, data.numAtts, data.varAtts, " " );MB_CHK_SET_ERR( rval, "Trouble getting attributes for variable " << data.varName ); 00418 } 00419 00420 return MB_SUCCESS; 00421 } 00422 00423 ErrorCode ReadNC::read_tag_values( const char*, const char*, const FileOptions&, std::vector< int >&, 00424 const SubsetList* ) 00425 { 00426 return MB_FAILURE; 00427 } 00428 00429 } // namespace moab