MOAB: Mesh Oriented datABase
(version 5.2.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 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