moab
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 #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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines