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