moab
NCHelper.cpp
Go to the documentation of this file.
00001 #include "NCHelper.hpp"
00002 #include "NCHelperEuler.hpp"
00003 #include "NCHelperFV.hpp"
00004 #include "NCHelperHOMME.hpp"
00005 #include "NCHelperMPAS.hpp"
00006 
00007 #include <sstream>
00008 
00009 #include "moab/ReadUtilIface.hpp"
00010 #include "MBTagConventions.hpp"
00011 
00012 #define ERRORR(rval, str) \
00013   if (MB_SUCCESS != rval) {_readNC->readMeshIface->report_error("%s", str); return rval;}
00014 
00015 #define ERRORS(err, str) \
00016   if (err) {_readNC->readMeshIface->report_error("%s", str); return MB_FAILURE;}
00017 
00018 namespace moab {
00019 
00020 NCHelper* NCHelper::get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
00021 {
00022   // Check if CF convention is being followed
00023   bool is_CF = false;
00024 
00025   std::map<std::string, ReadNC::AttData>& globalAtts = readNC->globalAtts;
00026   std::map<std::string, ReadNC::AttData>::iterator attIt = globalAtts.find("conventions");
00027   if (attIt == globalAtts.end())
00028     attIt = globalAtts.find("Conventions");
00029 
00030   if (attIt != globalAtts.end()) {
00031     unsigned int sz = attIt->second.attLen;
00032     std::string att_data;
00033     att_data.resize(sz + 1);
00034     att_data[sz] = '\000';
00035     int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
00036     if (0 == success && att_data.find("CF") != std::string::npos)
00037       is_CF = true;
00038   }
00039 
00040   if (is_CF) {
00041     if (NCHelperEuler::can_read_file(readNC, fileId))
00042       return new (std::nothrow) NCHelperEuler(readNC, fileId, opts, fileSet);
00043     else if (NCHelperFV::can_read_file(readNC, fileId))
00044       return new (std::nothrow) NCHelperFV(readNC, fileId, opts, fileSet);
00045     else if (NCHelperHOMME::can_read_file(readNC, fileId))
00046       return new (std::nothrow) NCHelperHOMME(readNC, fileId, opts, fileSet);
00047   }
00048   else {
00049     if (NCHelperMPAS::can_read_file(readNC))
00050       return new (std::nothrow) NCHelperMPAS(readNC, fileId, opts, fileSet);
00051     // For a HOMME connectivity file, there might be no CF convention
00052     else if (NCHelperHOMME::can_read_file(readNC, fileId))
00053       return new (std::nothrow) NCHelperHOMME(readNC, fileId, opts, fileSet);
00054   }
00055 
00056   // Unknown NetCDF grid (will fill this in later for POP, CICE and CLM)
00057   return NULL;
00058 }
00059 
00060 ErrorCode NCHelper::create_conventional_tags(const std::vector<int>& tstep_nums) {
00061   Interface*& mbImpl = _readNC->mbImpl;
00062   std::vector<std::string>& dimNames = _readNC->dimNames;
00063   std::vector<int>& dimLens = _readNC->dimLens;
00064   std::map<std::string, ReadNC::AttData>& globalAtts = _readNC->globalAtts;
00065   std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
00066   DebugOutput& dbgOut = _readNC->dbgOut;
00067   int& partMethod = _readNC->partMethod;
00068   ScdInterface* scdi = _readNC->scdi;
00069 
00070   ErrorCode rval;
00071   std::string tag_name;
00072 
00073   // <__NUM_DIMS>
00074   Tag numDimsTag = 0;
00075   tag_name = "__NUM_DIMS";
00076   int numDims = dimNames.size();
00077   rval = mbImpl->tag_get_handle(tag_name.c_str(), 1, MB_TYPE_INTEGER, numDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT);
00078   ERRORR(rval, "Trouble creating __NUM_DIMS tag.");
00079   rval = mbImpl->tag_set_data(numDimsTag, &_fileSet, 1, &numDims);
00080   ERRORR(rval, "Trouble setting data for __NUM_DIMS tag.");
00081   if (MB_SUCCESS == rval)
00082     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00083 
00084   // <__NUM_VARS>
00085   Tag numVarsTag = 0;
00086   tag_name = "__NUM_VARS";
00087   int numVars = varInfo.size();
00088   rval = mbImpl->tag_get_handle(tag_name.c_str(), 1, MB_TYPE_INTEGER, numVarsTag, MB_TAG_SPARSE | MB_TAG_CREAT);
00089   ERRORR(rval, "Trouble creating __NUM_VARS tag.");
00090   rval = mbImpl->tag_set_data(numVarsTag, &_fileSet, 1, &numVars);
00091   ERRORR(rval, "Trouble setting data for __NUM_VARS tag.");
00092   if (MB_SUCCESS == rval)
00093     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00094 
00095   // <__DIM_NAMES>
00096   Tag dimNamesTag = 0;
00097   tag_name = "__DIM_NAMES";
00098   std::string dimnames;
00099   unsigned int dimNamesSz = dimNames.size();
00100   for (unsigned int i = 0; i != dimNamesSz; i++) {
00101     dimnames.append(dimNames[i]);
00102     dimnames.push_back('\0');
00103   }
00104   int dimnamesSz = dimnames.size();
00105   rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, dimNamesTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00106   ERRORR(rval, "Trouble creating __DIM_NAMES tag.");
00107   const void* ptr = dimnames.c_str();
00108   rval = mbImpl->tag_set_by_ptr(dimNamesTag, &_fileSet, 1, &ptr, &dimnamesSz);
00109   ERRORR(rval, "Trouble setting data for __DIM_NAMES tag.");
00110   if (MB_SUCCESS == rval)
00111     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00112 
00113   // <__DIM_LENS>
00114   Tag dimLensTag = 0;
00115   tag_name = "__DIM_LENS";
00116   int dimLensSz = dimLens.size();
00117   rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_INTEGER, dimLensTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00118   ERRORR(rval, "Trouble creating __DIM_LENS tag.");
00119   ptr = &(dimLens[0]);
00120   rval = mbImpl->tag_set_by_ptr(dimLensTag, &_fileSet, 1, &ptr, &dimLensSz);
00121   ERRORR(rval, "Trouble setting data for __DIM_LENS tag.");
00122   if (MB_SUCCESS == rval)
00123     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00124 
00125   // <__VAR_NAMES>
00126   Tag varNamesTag = 0;
00127   tag_name = "__VAR_NAMES";
00128   std::string varnames;
00129   std::map<std::string, ReadNC::VarData>::iterator mapIter;
00130   for (mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter) {
00131     varnames.append(mapIter->first);
00132     varnames.push_back('\0');
00133   }
00134   int varnamesSz = varnames.size();
00135   rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, varNamesTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00136   ERRORR(rval, "Trouble creating __VAR_NAMES tag.");
00137   ptr = varnames.c_str();
00138   rval = mbImpl->tag_set_by_ptr(varNamesTag, &_fileSet, 1, &ptr, &varnamesSz);
00139   ERRORR(rval, "Trouble setting data for __VAR_NAMES tag.");
00140   if (MB_SUCCESS == rval)
00141     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00142 
00143   // __<dim_name>_LOC_MINMAX (for time)
00144   for (unsigned int i = 0; i != dimNamesSz; i++) {
00145     if (dimNames[i] == "time" || dimNames[i] == "Time" || dimNames[i] == "t") {
00146       std::stringstream ss_tag_name;
00147       ss_tag_name << "__" << dimNames[i] << "_LOC_MINMAX";
00148       tag_name = ss_tag_name.str();
00149       Tag tagh = 0;
00150       std::vector<int> val(2, 0);
00151       val[0] = 0;
00152       val[1] = nTimeSteps - 1;
00153       rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
00154       ERRORR(rval, "Trouble creating __<dim_name>_LOC_MINMAX tag.");
00155       rval = mbImpl->tag_set_data(tagh, &_fileSet, 1, &val[0]);
00156       ERRORR(rval, "Trouble setting data for __<dim_name>_LOC_MINMAX tag.");
00157       if (MB_SUCCESS == rval)
00158         dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00159     }
00160   }
00161 
00162   // __<dim_name>_LOC_VALS (for time)
00163   for (unsigned int i = 0; i != dimNamesSz; i++) {
00164     if (dimNames[i] == "time" || dimNames[i] == "Time" || dimNames[i] == "t") {
00165       std::vector<int> val;
00166       if (!tstep_nums.empty())
00167         val = tstep_nums;
00168       else {
00169         val.resize(tVals.size());
00170         for (unsigned int j = 0; j != tVals.size(); j++)
00171           val[j] = j;
00172       }
00173       Tag tagh = 0;
00174       std::stringstream ss_tag_name;
00175       ss_tag_name << "__" << dimNames[i] << "_LOC_VALS";
00176       tag_name = ss_tag_name.str();
00177       rval = mbImpl->tag_get_handle(tag_name.c_str(), val.size(), MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
00178       ERRORR(rval, "Trouble creating __<dim_name>_LOC_VALS tag.");
00179       rval = mbImpl->tag_set_data(tagh, &_fileSet, 1, &val[0]);
00180       ERRORR(rval, "Trouble setting data for __<dim_name>_LOC_VALS tag.");
00181       if (MB_SUCCESS == rval)
00182         dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00183     }
00184   }
00185 
00186   // __<var_name>_DIMS
00187   for (mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter) {
00188     Tag varNamesDimsTag = 0;
00189     std::stringstream ss_tag_name;
00190     ss_tag_name << "__" << mapIter->first << "_DIMS";
00191     tag_name = ss_tag_name.str();
00192     unsigned int varDimSz = varInfo[mapIter->first].varDims.size();
00193     if (varDimSz == 0)
00194       continue;
00195     varInfo[mapIter->first].varTags.resize(varDimSz, 0);
00196     for (unsigned int i = 0; i != varDimSz; i++) {
00197       Tag tmptag = 0;
00198       std::string tmptagname = dimNames[varInfo[mapIter->first].varDims[i]];
00199       mbImpl->tag_get_handle(tmptagname.c_str(), 0, MB_TYPE_OPAQUE, tmptag, MB_TAG_ANY);
00200       varInfo[mapIter->first].varTags[i] = tmptag;
00201     }
00202     rval = mbImpl->tag_get_handle(tag_name.c_str(), varDimSz, MB_TYPE_HANDLE, varNamesDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT);
00203     ERRORR(rval, "Trouble creating __<var_name>_DIMS tag.");
00204     rval = mbImpl->tag_set_data(varNamesDimsTag, &_fileSet, 1, &(varInfo[mapIter->first].varTags[0]));
00205     ERRORR(rval, "Trouble setting data for __<var_name>_DIMS tag.");
00206     if (MB_SUCCESS == rval)
00207       dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00208   }
00209 
00210   // <PARTITION_METHOD>
00211   Tag part_tag = scdi->part_method_tag();
00212   if (!part_tag)
00213     ERRORR(MB_FAILURE, "Trouble getting partition method tag.");
00214   rval = mbImpl->tag_set_data(part_tag, &_fileSet, 1, &partMethod);
00215   ERRORR(rval, "Trouble setting data for PARTITION_METHOD tag.");
00216   if (MB_SUCCESS == rval)
00217     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00218 
00219   // <__GLOBAL_ATTRIBS>
00220   tag_name = "__GLOBAL_ATTRIBS";
00221   Tag globalAttTag = 0;
00222   rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, globalAttTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00223   ERRORR(rval, "Trouble creating __GLOBAL_ATTRIBS tag.");
00224   std::string gattVal;
00225   std::vector<int> gattLen;
00226   rval = create_attrib_string(globalAtts, gattVal, gattLen);
00227   ERRORR(rval, "Trouble creating attribute strings.");
00228   const void* gattptr = gattVal.c_str();
00229   int globalAttSz = gattVal.size();
00230   rval = mbImpl->tag_set_by_ptr(globalAttTag, &_fileSet, 1, &gattptr, &globalAttSz);
00231   ERRORR(rval, "Trouble setting data for __GLOBAL_ATTRIBS tag.");
00232   if (MB_SUCCESS == rval)
00233     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00234 
00235   // <__GLOBAL_ATTRIBS_LEN>
00236   tag_name = "__GLOBAL_ATTRIBS_LEN";
00237   Tag globalAttLenTag = 0;
00238   if (gattLen.size() == 0)
00239     gattLen.push_back(0);
00240   rval = mbImpl->tag_get_handle(tag_name.c_str(), gattLen.size(), MB_TYPE_INTEGER, globalAttLenTag, MB_TAG_SPARSE | MB_TAG_CREAT);
00241   ERRORR(rval, "Trouble creating __GLOBAL_ATTRIBS_LEN tag.");
00242   rval = mbImpl->tag_set_data(globalAttLenTag, &_fileSet, 1, &gattLen[0]);
00243   ERRORR(rval, "Trouble setting data for __GLOBAL_ATTRIBS_LEN tag.");
00244   if (MB_SUCCESS == rval)
00245     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00246 
00247   // __<var_name>_ATTRIBS and __<var_name>_ATTRIBS_LEN
00248   for (mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter) {
00249     std::stringstream ssTagName;
00250     ssTagName << "__" << mapIter->first << "_ATTRIBS";
00251     tag_name = ssTagName.str();
00252     Tag varAttTag = 0;
00253     rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, varAttTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00254     ERRORR(rval, "Trouble creating __<var_name>_ATTRIBS tag.");
00255     std::string varAttVal;
00256     std::vector<int> varAttLen;
00257     rval = create_attrib_string(mapIter->second.varAtts, varAttVal, varAttLen);
00258     ERRORR(rval, "Trouble creating attribute strings.");
00259     const void* varAttPtr = varAttVal.c_str();
00260     int varAttSz = varAttVal.size();
00261     rval = mbImpl->tag_set_by_ptr(varAttTag, &_fileSet, 1, &varAttPtr, &varAttSz);
00262     ERRORR(rval, "Trouble setting data for __<var_name>_ATTRIBS tag.");
00263     if (MB_SUCCESS == rval)
00264       dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00265     if (varAttLen.size() == 0)
00266       varAttLen.push_back(0);
00267     ssTagName << "_LEN";
00268     tag_name = ssTagName.str();
00269     Tag varAttLenTag = 0;
00270     rval = mbImpl->tag_get_handle(tag_name.c_str(), varAttLen.size(), MB_TYPE_INTEGER, varAttLenTag, MB_TAG_SPARSE | MB_TAG_CREAT);
00271     ERRORR(rval, "Trouble creating __<var_name>_ATTRIBS_LEN tag.");
00272     rval = mbImpl->tag_set_data(varAttLenTag, &_fileSet, 1, &varAttLen[0]);
00273     ERRORR(rval, "Trouble setting data for __<var_name>_ATTRIBS_LEN tag.");
00274     if (MB_SUCCESS == rval)
00275       dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00276   }
00277 
00278   // <__VAR_NAMES_LOCATIONS>
00279   tag_name = "__VAR_NAMES_LOCATIONS";
00280   Tag varNamesLocsTag = 0;
00281   std::vector<int> varNamesLocs(varInfo.size());
00282   rval = mbImpl->tag_get_handle(tag_name.c_str(), varNamesLocs.size(), MB_TYPE_INTEGER, varNamesLocsTag, MB_TAG_CREAT
00283       | MB_TAG_SPARSE);
00284   ERRORR(rval, "Trouble creating __VAR_NAMES_LOCATIONS tag.");
00285   for (mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter) {
00286     varNamesLocs[std::distance(varInfo.begin(), mapIter)] = mapIter->second.entLoc;
00287   }
00288   rval = mbImpl->tag_set_data(varNamesLocsTag, &_fileSet, 1, &varNamesLocs[0]);
00289   ERRORR(rval, "Trouble setting data for __VAR_NAMES_LOCATIONS tag.");
00290   if (MB_SUCCESS == rval)
00291     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00292 
00293   // <__MESH_TYPE>
00294   Tag meshTypeTag = 0;
00295   tag_name = "__MESH_TYPE";
00296   std::string meshTypeName = get_mesh_type_name();
00297 
00298   rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, meshTypeTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00299   ERRORR(rval, "Trouble creating __MESH_TYPE tag.");
00300   ptr = meshTypeName.c_str();
00301   int leng = meshTypeName.size();
00302   rval = mbImpl->tag_set_by_ptr(meshTypeTag, &_fileSet, 1, &ptr, &leng);
00303   ERRORR(rval, "Trouble setting data for __MESH_TYPE tag.");
00304   if (MB_SUCCESS == rval)
00305     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
00306 
00307   return MB_SUCCESS;
00308 }
00309 
00310 ErrorCode NCHelper::read_variable_setup(std::vector<std::string>& var_names, std::vector<int>& tstep_nums,
00311                                         std::vector<ReadNC::VarData>& vdatas, std::vector<ReadNC::VarData>& vsetdatas)
00312 {
00313   std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
00314   std::map<std::string, ReadNC::VarData>::iterator mit;
00315 
00316   // If empty read them all (except ignored variables)
00317   if (var_names.empty()) {
00318     for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) {
00319       ReadNC::VarData vd = (*mit).second;
00320 
00321       // If read all variables at once, skip ignored ones
00322       // Upon creation of dummy variables, tag values have already been set
00323       if (ignoredVarNames.find(vd.varName) != ignoredVarNames.end() ||
00324           dummyVarNames.find(vd.varName) != dummyVarNames.end())
00325          continue;
00326 
00327       if (vd.entLoc == ReadNC::ENTLOCSET)
00328         vsetdatas.push_back(vd);
00329       else
00330         vdatas.push_back(vd);
00331     }
00332   }
00333   else {
00334     // Read specified variables (might include ignored ones)
00335     for (unsigned int i = 0; i < var_names.size(); i++) {
00336       mit = varInfo.find(var_names[i]);
00337       if (mit != varInfo.end()) {
00338         ReadNC::VarData vd = (*mit).second;
00339 
00340         // Upon creation of dummy variables, tag values have already been set
00341         if (dummyVarNames.find(vd.varName) != dummyVarNames.end())
00342            continue;
00343 
00344         if (vd.entLoc == ReadNC::ENTLOCSET)
00345           vsetdatas.push_back(vd);
00346         else
00347           vdatas.push_back(vd);
00348       }
00349       else {
00350         ERRORR(MB_FAILURE, "Couldn't find variable.");
00351       }
00352     }
00353   }
00354 
00355   if (tstep_nums.empty() && nTimeSteps > 0) {
00356     // No timesteps input, get them all
00357     for (int i = 0; i < nTimeSteps; i++)
00358       tstep_nums.push_back(i);
00359   }
00360 
00361   if (!tstep_nums.empty()) {
00362     for (unsigned int i = 0; i < vdatas.size(); i++) {
00363       vdatas[i].varTags.resize(tstep_nums.size(), 0);
00364       vdatas[i].varDatas.resize(tstep_nums.size());
00365       vdatas[i].has_tsteps = true;
00366     }
00367 
00368     for (unsigned int i = 0; i < vsetdatas.size(); i++) {
00369       if ((std::find(vsetdatas[i].varDims.begin(), vsetdatas[i].varDims.end(), tDim) != vsetdatas[i].varDims.end())
00370           && (vsetdatas[i].varDims.size() > 1)) {
00371         // Set variables with timesteps: time is the first dimension, followed
00372         // by other dimensions, e.g. xtime(Time, StrLen)
00373         vsetdatas[i].varTags.resize(tstep_nums.size(), 0);
00374         vsetdatas[i].varDatas.resize(tstep_nums.size());
00375         vsetdatas[i].has_tsteps = true;
00376       }
00377       else {
00378         // Set variables without timesteps: no time dimension, or time is the only
00379         // dimension, e.g. lev(lev), xtime(Time)
00380         vsetdatas[i].varTags.resize(1, 0);
00381         vsetdatas[i].varDatas.resize(1);
00382         vsetdatas[i].has_tsteps = false;
00383       }
00384     }
00385   }
00386 
00387   return MB_SUCCESS;
00388 }
00389 
00390 ErrorCode NCHelper::read_variable_to_set(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums)
00391 {
00392   Interface*& mbImpl = _readNC->mbImpl;
00393   DebugOutput& dbgOut = _readNC->dbgOut;
00394 
00395   ErrorCode rval = read_variable_to_set_allocate(vdatas, tstep_nums);
00396   ERRORR(rval, "Trouble allocating read variables to set.");
00397 
00398   // Finally, read into that space
00399   int success;
00400   for (unsigned int i = 0; i < vdatas.size(); i++) {
00401     // Note, for set variables without timesteps, loop one time and then break
00402     for (unsigned int t = 0; t < tstep_nums.size(); t++) {
00403       void* data = vdatas[i].varDatas[t];
00404 
00405       // Set variables with timesteps, e.g. xtime(Time, StrLen)
00406       if (vdatas[i].has_tsteps) {
00407         // Set readStart for each timestep along time dimension
00408         vdatas[i].readStarts[0] = tstep_nums[t];
00409       }
00410 
00411       switch (vdatas[i].varDataType) {
00412         case NC_BYTE:
00413         case NC_CHAR:
00414           success = NCFUNCAG(_vara_text)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
00415                                         &vdatas[i].readCounts[0], (char*) data);
00416           ERRORS(success, "Failed to read char data.");
00417           break;
00418         case NC_DOUBLE:
00419           success = NCFUNCAG(_vara_double)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
00420                                         &vdatas[i].readCounts[0], (double*) data);
00421           ERRORS(success, "Failed to read double data.");
00422           break;
00423         case NC_FLOAT:
00424           success = NCFUNCAG(_vara_float)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
00425                                         &vdatas[i].readCounts[0], (float*) data);
00426           ERRORS(success, "Failed to read float data.");
00427           break;
00428         case NC_INT:
00429           success = NCFUNCAG(_vara_int)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
00430                                         &vdatas[i].readCounts[0], (int*) data);
00431           ERRORS(success, "Failed to read int data.");
00432           break;
00433         case NC_SHORT:
00434           success = NCFUNCAG(_vara_short)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
00435                                         &vdatas[i].readCounts[0], (short*) data);
00436           ERRORS(success, "Failed to read short data.");
00437           break;
00438         default:
00439           success = 1;
00440       }
00441 
00442       if (success)
00443         ERRORR(MB_FAILURE, "Trouble reading variable.");
00444 
00445       dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
00446       rval = convert_variable(vdatas[i], t);
00447       ERRORR(rval, "Failed to convert variable.");
00448 
00449       dbgOut.tprintf(2, "Setting data for variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
00450       rval = mbImpl->tag_set_by_ptr(vdatas[i].varTags[t], &_fileSet, 1, &data, &vdatas[i].sz);
00451       ERRORR(rval, "Failed to set data for variable.");
00452 
00453       // Memory pointed by pointer data can be deleted, as tag_set_by_ptr() has already copied the tag values
00454       switch (vdatas[i].varDataType) {
00455         case NC_BYTE:
00456         case NC_CHAR:
00457           delete[] (char*) data;
00458           break;
00459         case NC_DOUBLE:
00460         case NC_FLOAT:
00461           delete[] (double*) data;
00462           break;
00463         case NC_INT:
00464         case NC_SHORT:
00465           delete[] (int*) data;
00466           break;
00467         default:
00468           break;
00469       }
00470       vdatas[i].varDatas[t] = NULL;
00471 
00472       // Loop continues only for set variables with timesteps, e.g. xtime(Time, StrLen)
00473       if (!vdatas[i].has_tsteps)
00474         break;
00475     }
00476   }
00477 
00478   // Debug output, if requested
00479   if (1 == dbgOut.get_verbosity()) {
00480     dbgOut.printf(1, "Read variables: %s", vdatas.begin()->varName.c_str());
00481     for (unsigned int i = 1; i < vdatas.size(); i++)
00482       dbgOut.printf(1, ", %s ", vdatas[i].varName.c_str());
00483     dbgOut.tprintf(1, "\n");
00484   }
00485 
00486   return rval;
00487 }
00488 
00489 ErrorCode NCHelper::convert_variable(ReadNC::VarData& var_data, int tstep_num)
00490 {
00491   DebugOutput& dbgOut = _readNC->dbgOut;
00492 
00493   // Get ptr to tag space
00494   void* data = var_data.varDatas[tstep_num];
00495 
00496   // Get variable size
00497   std::size_t sz = var_data.sz;
00498   assert(sz > 0);
00499 
00500   // Finally, read into that space
00501   int success = 0;
00502   int* idata;
00503   double* ddata;
00504   float* fdata;
00505   short* sdata;
00506 
00507   switch (var_data.varDataType) {
00508     case NC_FLOAT:
00509       ddata = (double*) var_data.varDatas[tstep_num];
00510       fdata = (float*) var_data.varDatas[tstep_num];
00511       // Convert in-place
00512       for (int i = sz - 1; i >= 0; i--)
00513         ddata[i] = fdata[i];
00514       break;
00515     case NC_SHORT:
00516       idata = (int*) var_data.varDatas[tstep_num];
00517       sdata = (short*) var_data.varDatas[tstep_num];
00518       // Convert in-place
00519       for (int i = sz - 1; i >= 0; i--)
00520         idata[i] = sdata[i];
00521       break;
00522     default:
00523       success = 1;
00524   }
00525 
00526   if (2 <= dbgOut.get_verbosity() && !success) {
00527     double dmin, dmax;
00528     int imin, imax;
00529     switch (var_data.varDataType) {
00530       case NC_DOUBLE:
00531       case NC_FLOAT:
00532         ddata = (double*) data;
00533         if (sz == 0)
00534           break;
00535 
00536         dmin = dmax = ddata[0];
00537         for (unsigned int i = 1; i < sz; i++) {
00538           if (ddata[i] < dmin)
00539             dmin = ddata[i];
00540           if (ddata[i] > dmax)
00541             dmax = ddata[i];
00542         }
00543         dbgOut.tprintf(2, "Variable %s (double): min = %f, max = %f\n", var_data.varName.c_str(), dmin, dmax);
00544         break;
00545       case NC_INT:
00546       case NC_SHORT:
00547         idata = (int*) data;
00548         if (sz == 0)
00549           break;
00550 
00551         imin = imax = idata[0];
00552         for (unsigned int i = 1; i < sz; i++) {
00553           if (idata[i] < imin)
00554             imin = idata[i];
00555           if (idata[i] > imax)
00556             imax = idata[i];
00557         }
00558         dbgOut.tprintf(2, "Variable %s (int): min = %d, max = %d\n", var_data.varName.c_str(), imin, imax);
00559         break;
00560       case NC_NAT:
00561       case NC_BYTE:
00562       case NC_CHAR:
00563         break;
00564       default: // Default case added to remove compiler warnings
00565         success = 1;
00566     }
00567   }
00568 
00569   return MB_SUCCESS;
00570 }
00571 
00572 ErrorCode NCHelper::read_coordinate(const char* var_name, int lmin, int lmax, std::vector<double>& cvals)
00573 {
00574   std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
00575   std::map<std::string, ReadNC::VarData>::iterator vmit = varInfo.find(var_name);
00576   if (varInfo.end() == vmit)
00577     return MB_FAILURE;
00578 
00579   assert(lmin >= 0 && lmax >= lmin);
00580   NCDF_SIZE tstart = lmin;
00581   NCDF_SIZE tcount = lmax - lmin + 1;
00582   NCDF_DIFF dum_stride = 1;
00583   int fail;
00584 
00585   // Check size
00586   if ((std::size_t)tcount != cvals.size())
00587     cvals.resize(tcount);
00588 
00589   // Check to make sure it's a float or double
00590   if (NC_DOUBLE == (*vmit).second.varDataType) {
00591     fail = NCFUNCAG(_vars_double)(_fileId, (*vmit).second.varId, &tstart, &tcount, &dum_stride, &cvals[0]);
00592     if (fail)
00593       ERRORS(MB_FAILURE, "Failed to get coordinate values.");
00594   }
00595   else if (NC_FLOAT == (*vmit).second.varDataType) {
00596     std::vector<float> tcvals(tcount);
00597     fail = NCFUNCAG(_vars_float)(_fileId, (*vmit).second.varId, &tstart, &tcount, &dum_stride, &tcvals[0]);
00598     if (fail)
00599       ERRORS(MB_FAILURE, "Failed to get coordinate values.");
00600     std::copy(tcvals.begin(), tcvals.end(), cvals.begin());
00601   }
00602   else {
00603     ERRORR(MB_FAILURE, "Wrong data type for coordinate variable.");
00604   }
00605 
00606   return MB_SUCCESS;
00607 }
00608 
00609 ErrorCode NCHelper::get_tag_to_set(ReadNC::VarData& var_data, int tstep_num, Tag& tagh)
00610 {
00611   Interface*& mbImpl = _readNC->mbImpl;
00612   DebugOutput& dbgOut = _readNC->dbgOut;
00613 
00614   std::ostringstream tag_name;
00615   if (var_data.has_tsteps)
00616     tag_name << var_data.varName << tstep_num;
00617   else
00618     tag_name << var_data.varName;
00619 
00620   ErrorCode rval = MB_SUCCESS;
00621   tagh = 0;
00622   switch (var_data.varDataType) {
00623     case NC_BYTE:
00624     case NC_CHAR:
00625       rval = mbImpl->tag_get_handle(tag_name.str().c_str(), 0, MB_TYPE_OPAQUE, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00626       break;
00627     case NC_DOUBLE:
00628     case NC_FLOAT:
00629       rval = mbImpl->tag_get_handle(tag_name.str().c_str(), 0, MB_TYPE_DOUBLE, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00630       break;
00631     case NC_INT:
00632     case NC_SHORT:
00633       rval = mbImpl->tag_get_handle(tag_name.str().c_str(), 0, MB_TYPE_INTEGER, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00634       break;
00635     default:
00636       std::cerr << "Unrecognized data type for tag " << tag_name << std::endl;
00637       rval = MB_FAILURE;
00638   }
00639 
00640   if (MB_SUCCESS == rval)
00641     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.str().c_str());
00642 
00643   return rval;
00644 }
00645 
00646 ErrorCode NCHelper::get_tag_to_nonset(ReadNC::VarData& var_data, int tstep_num, Tag& tagh, int num_lev)
00647 {
00648   Interface*& mbImpl = _readNC->mbImpl;
00649   DebugOutput& dbgOut = _readNC->dbgOut;
00650 
00651   std::ostringstream tag_name;
00652   tag_name << var_data.varName << tstep_num;
00653 
00654   ErrorCode rval = MB_SUCCESS;
00655   tagh = 0;
00656   switch (var_data.varDataType) {
00657     case NC_BYTE:
00658     case NC_CHAR:
00659       rval = mbImpl->tag_get_handle(tag_name.str().c_str(), num_lev, MB_TYPE_OPAQUE, tagh, MB_TAG_DENSE | MB_TAG_CREAT);
00660       break;
00661     case NC_DOUBLE:
00662     case NC_FLOAT:
00663       rval = mbImpl->tag_get_handle(tag_name.str().c_str(), num_lev, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT);
00664       break;
00665     case NC_INT:
00666     case NC_SHORT:
00667       rval = mbImpl->tag_get_handle(tag_name.str().c_str(), num_lev, MB_TYPE_INTEGER, tagh, MB_TAG_DENSE | MB_TAG_CREAT);
00668       break;
00669     default:
00670       std::cerr << "Unrecognized data type for tag " << tag_name.str() << std::endl;
00671       rval = MB_FAILURE;
00672   }
00673 
00674   if (MB_SUCCESS == rval)
00675     dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.str().c_str());
00676 
00677   return rval;
00678 }
00679 
00680 ErrorCode NCHelper::create_attrib_string(const std::map<std::string, ReadNC::AttData>& attMap, std::string& attVal, std::vector<int>& attLen)
00681 {
00682   int success;
00683   std::stringstream ssAtt;
00684   unsigned int sz = 0;
00685   std::map<std::string, ReadNC::AttData>::const_iterator attIt = attMap.begin();
00686   for (; attIt != attMap.end(); ++attIt) {
00687     ssAtt << attIt->second.attName;
00688     ssAtt << '\0';
00689     void* attData = NULL;
00690     switch (attIt->second.attDataType) {
00691       case NC_BYTE:
00692       case NC_CHAR:
00693         sz = attIt->second.attLen;
00694         attData = (char *) malloc(sz);
00695         success = NCFUNC(get_att_text)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (char*) attData);
00696         ERRORS(success, "Failed to read attribute char data.");
00697         ssAtt << "char;";
00698         break;
00699       case NC_DOUBLE:
00700         sz = attIt->second.attLen * sizeof(double);
00701         attData = (double *) malloc(sz);
00702         success = NCFUNC(get_att_double)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (double*) attData);
00703         ERRORS(success, "Failed to read attribute double data.");
00704         ssAtt << "double;";
00705         break;
00706       case NC_FLOAT:
00707         sz = attIt->second.attLen * sizeof(float);
00708         attData = (float *) malloc(sz);
00709         success = NCFUNC(get_att_float)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (float*) attData);
00710         ERRORS(success, "Failed to read attribute float data.");
00711         ssAtt << "float;";
00712         break;
00713       case NC_INT:
00714         sz = attIt->second.attLen * sizeof(int);
00715         attData = (int *) malloc(sz);
00716         success = NCFUNC(get_att_int)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (int*) attData);
00717         ERRORS(success, "Failed to read attribute int data.");
00718         ssAtt << "int;";
00719         break;
00720       case NC_SHORT:
00721         sz = attIt->second.attLen * sizeof(short);
00722         attData = (short *) malloc(sz);
00723         success = NCFUNC(get_att_short)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (short*) attData);
00724         ERRORS(success, "Failed to read attribute short data.");
00725         ssAtt << "short;";
00726         break;
00727       default:
00728         success = 1;
00729     }
00730     char* tmpc = (char *) attData;
00731     for (unsigned int counter = 0; counter != sz; ++counter)
00732       ssAtt << tmpc[counter];
00733     free(attData);
00734     ssAtt << ';';
00735     attLen.push_back(ssAtt.str().size() - 1);
00736   }
00737   attVal = ssAtt.str();
00738 
00739   return MB_SUCCESS;
00740 }
00741 
00742 ErrorCode NCHelper::create_dummy_variables()
00743 {
00744   Interface*& mbImpl = _readNC->mbImpl;
00745   std::vector<std::string>& dimNames = _readNC->dimNames;
00746   std::vector<int>& dimLens = _readNC->dimLens;
00747   std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
00748   DebugOutput& dbgOut = _readNC->dbgOut;
00749 
00750   // Hack: look at all dimensions, and see if we have one that does not appear in the list of varInfo names
00751   // Right now, candidates are from unstructured meshes, such as ncol (HOMME) and nCells (MPAS)
00752   // For each of them, create a dummy variable with a sparse tag to store the dimension length
00753   for (unsigned int i = 0; i < dimNames.size(); i++) {
00754     // If there is a variable with this dimension name, skip
00755     if (varInfo.find(dimNames[i]) != varInfo.end())
00756       continue;
00757 
00758     // Create a dummy variable
00759     int sizeTotalVar = varInfo.size();
00760     std::string var_name(dimNames[i]);
00761     ReadNC::VarData& data = varInfo[var_name];
00762     data.varName = std::string(var_name);
00763     data.varId = sizeTotalVar;
00764     data.varTags.resize(1, 0);
00765     data.varDataType = NC_INT;
00766     data.varDims.resize(1);
00767     data.varDims[0] = (int)i;
00768     data.numAtts = 0;
00769     data.entLoc = ReadNC::ENTLOCSET;
00770     dummyVarNames.insert(dimNames[i]);
00771     dbgOut.tprintf(2, "Dummy variable created for dimension %s\n", dimNames[i].c_str());
00772 
00773     // Create a corresponding sparse tag
00774     Tag tagh;
00775     ErrorCode rval = mbImpl->tag_get_handle(dimNames[i].c_str(), 0, MB_TYPE_INTEGER, tagh,
00776                                             MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
00777     ERRORR(rval, "Failed to create tag for a dummy dimension variable.");
00778 
00779     // Tag value is the dimension length
00780     const void* ptr = &dimLens[i];
00781     // Tag size is 1
00782     int size = 1;
00783     rval = mbImpl->tag_set_by_ptr(tagh, &_fileSet, 1, &ptr, &size);
00784     ERRORR(rval, "Failed to set data for dimension tag.");
00785 
00786     dbgOut.tprintf(2, "Sparse tag created for dimension %s\n", dimNames[i].c_str());
00787   }
00788 
00789   return MB_SUCCESS;
00790 }
00791 
00792 ErrorCode NCHelper::read_variable_to_set_allocate(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums)
00793 {
00794   std::vector<int>& dimLens = _readNC->dimLens;
00795   DebugOutput& dbgOut = _readNC->dbgOut;
00796 
00797   ErrorCode rval = MB_SUCCESS;
00798 
00799   for (unsigned int i = 0; i < vdatas.size(); i++) {
00800     // Set up readStarts and readCounts
00801     if (vdatas[i].has_tsteps) {
00802       // First: time
00803       vdatas[i].readStarts.push_back(0); // This value is timestep dependent, will be set later
00804       vdatas[i].readCounts.push_back(1);
00805 
00806       // Next: other dimensions
00807       for (unsigned int idx = 1; idx != vdatas[i].varDims.size(); idx++){
00808         vdatas[i].readStarts.push_back(0);
00809         vdatas[i].readCounts.push_back(dimLens[vdatas[i].varDims[idx]]);
00810       }
00811     }
00812     else {
00813       if (vdatas[i].varDims.empty()) {
00814         // Scalar variable
00815         vdatas[i].readStarts.push_back(0);
00816         vdatas[i].readCounts.push_back(1);
00817       }
00818       else {
00819         for (unsigned int idx = 0; idx != vdatas[i].varDims.size(); idx++){
00820           vdatas[i].readStarts.push_back(0);
00821           vdatas[i].readCounts.push_back(dimLens[vdatas[i].varDims[idx]]);
00822         }
00823       }
00824     }
00825 
00826     // Get variable size
00827     vdatas[i].sz = 1;
00828     for (std::size_t idx = 0; idx != vdatas[i].readCounts.size(); idx++)
00829       vdatas[i].sz *= vdatas[i].readCounts[idx];
00830 
00831     // Note, for set variables without timesteps, loop one time and then break
00832     for (unsigned int t = 0; t < tstep_nums.size(); t++) {
00833       dbgOut.tprintf(2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
00834 
00835       if (tstep_nums[t] >= dimLens[tDim]) {
00836         ERRORR(MB_INDEX_OUT_OF_RANGE, "Wrong value for a timestep number.");
00837       }
00838 
00839       // Get the tag to read into
00840       if (!vdatas[i].varTags[t]) {
00841         rval = get_tag_to_set(vdatas[i], tstep_nums[t], vdatas[i].varTags[t]);
00842         ERRORR(rval, "Trouble getting tag for a set variable.");
00843       }
00844 
00845       switch (vdatas[i].varDataType) {
00846         case NC_BYTE:
00847         case NC_CHAR:
00848           vdatas[i].varDatas[t] = new char[vdatas[i].sz];
00849           break;
00850         case NC_DOUBLE:
00851         case NC_FLOAT:
00852           vdatas[i].varDatas[t] = new double[vdatas[i].sz];
00853           break;
00854         case NC_INT:
00855         case NC_SHORT:
00856           vdatas[i].varDatas[t] = new int[vdatas[i].sz];
00857           break;
00858         default:
00859           std::cerr << "Unrecognized data type for set variable tag values" << std::endl;
00860           rval = MB_FAILURE;
00861       }
00862 
00863       // Loop continues only for set variables with timesteps, e.g. xtime(Time, StrLen)
00864       if (!vdatas[i].has_tsteps)
00865         break;
00866     }
00867   }
00868 
00869   return rval;
00870 }
00871 
00872 ErrorCode ScdNCHelper::check_existing_mesh() {
00873   Interface*& mbImpl = _readNC->mbImpl;
00874 
00875   // Get the number of vertices
00876   int num_verts;
00877   ErrorCode rval = mbImpl->get_number_entities_by_dimension(_fileSet, 0, num_verts);
00878   ERRORR(rval, "Trouble getting number of vertices.");
00879 
00880   /*
00881   // Check against parameters
00882   // When ghosting is used, this check might fail (to be updated later)
00883   if (num_verts > 0) {
00884     int expected_verts = (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1) * (-1 == lDims[2] ? 1 : lDims[5] - lDims[2] + 1);
00885     if (num_verts != expected_verts) {
00886       ERRORR(MB_FAILURE, "Number of vertices doesn't match.");
00887     }
00888   }
00889   */
00890 
00891   // Check the number of elements too
00892   int num_elems;
00893   rval = mbImpl->get_number_entities_by_dimension(_fileSet, (-1 == lCDims[2] ? 2 : 3), num_elems);
00894   ERRORR(rval, "Trouble getting number of elements.");
00895 
00896   /*
00897   // Check against parameters
00898   // When ghosting is used, this check might fail (to be updated later)
00899   if (num_elems > 0) {
00900     int expected_elems = (lCDims[3] - lCDims[0] + 1) * (lCDims[4] - lCDims[1] + 1) * (-1 == lCDims[2] ? 1 : (lCDims[5] - lCDims[2] + 1));
00901     if (num_elems != expected_elems) {
00902       ERRORR(MB_FAILURE, "Number of elements doesn't match.");
00903     }
00904   }
00905   */
00906 
00907   return MB_SUCCESS;
00908 }
00909 
00910 ErrorCode ScdNCHelper::create_mesh(Range& faces)
00911 {
00912   Interface*& mbImpl = _readNC->mbImpl;
00913   Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
00914   const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
00915   DebugOutput& dbgOut = _readNC->dbgOut;
00916   ScdInterface* scdi = _readNC->scdi;
00917   ScdParData& parData = _readNC->parData;
00918 
00919   Range tmp_range;
00920   ScdBox* scd_box;
00921 
00922   ErrorCode rval = scdi->construct_box(HomCoord(lDims[0], lDims[1], lDims[2], 1), HomCoord(lDims[3], lDims[4], lDims[5], 1), 
00923                                        NULL, 0, scd_box, locallyPeriodic, &parData, true);
00924   ERRORR(rval, "Trouble creating scd vertex sequence.");
00925 
00926   // Add verts to tmp_range first, so we can duplicate global ids in vertex ids
00927   tmp_range.insert(scd_box->start_vertex(), scd_box->start_vertex() + scd_box->num_vertices() - 1);
00928 
00929   if (mpFileIdTag) {
00930     int count;
00931     void* data;
00932     rval = mbImpl->tag_iterate(*mpFileIdTag, tmp_range.begin(), tmp_range.end(), count, data);
00933     ERRORR(rval, "Failed to get tag iterator on file id tag.");
00934     assert(count == scd_box->num_vertices());
00935     int* fid_data = (int*) data;
00936     rval = mbImpl->tag_iterate(mGlobalIdTag, tmp_range.begin(), tmp_range.end(), count, data);
00937     ERRORR(rval, "Failed to get tag iterator on global id tag.");
00938     assert(count == scd_box->num_vertices());
00939     int* gid_data = (int*) data;
00940     for (int i = 0; i < count; i++)
00941       fid_data[i] = gid_data[i];
00942   }
00943 
00944   // Then add box set and elements to the range, then to the file set
00945   tmp_range.insert(scd_box->start_element(), scd_box->start_element() + scd_box->num_elements() - 1);
00946   tmp_range.insert(scd_box->box_set());
00947   rval = mbImpl->add_entities(_fileSet, tmp_range);
00948   ERRORR(rval, "Couldn't add new vertices to file set.");
00949 
00950   dbgOut.tprintf(1, "scdbox %d quads, %d vertices\n", scd_box->num_elements(), scd_box->num_vertices());
00951 
00952   // Set the vertex coordinates
00953   double *xc, *yc, *zc;
00954   rval = scd_box->get_coordinate_arrays(xc, yc, zc);
00955   ERRORR(rval, "Couldn't get vertex coordinate arrays.");
00956 
00957   int i, j, k, il, jl, kl;
00958   int dil = lDims[3] - lDims[0] + 1;
00959   int djl = lDims[4] - lDims[1] + 1;
00960   assert(dil == (int)ilVals.size() && djl == (int)jlVals.size() &&
00961       (-1 == lDims[2] || lDims[5] - lDims[2] + 1 == (int)levVals.size()));
00962 //#define INDEX(i, j, k) ()
00963   for (kl = lDims[2]; kl <= lDims[5]; kl++) {
00964     k = kl - lDims[2];
00965     for (jl = lDims[1]; jl <= lDims[4]; jl++) {
00966       j = jl - lDims[1];
00967       for (il = lDims[0]; il <= lDims[3]; il++) {
00968         i = il - lDims[0];
00969         unsigned int pos = i + j * dil + k * dil * djl;
00970         xc[pos] = ilVals[i];
00971         yc[pos] = jlVals[j];
00972         zc[pos] = (-1 == lDims[2] ? 0.0 : levVals[k]);
00973       }
00974     }
00975   }
00976 //#undef INDEX
00977 
00978 #ifndef NDEBUG
00979   int num_verts = (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1) * (-1 == lDims[2] ? 1 : lDims[5] - lDims[2] + 1);
00980   std::vector<int> gids(num_verts);
00981   Range verts(scd_box->start_vertex(), scd_box->start_vertex() + scd_box->num_vertices() - 1);
00982   rval = mbImpl->tag_get_data(mGlobalIdTag, verts, &gids[0]);
00983   ERRORR(rval, "Trouble getting gid values.");
00984   int vmin = *(std::min_element(gids.begin(), gids.end())), vmax = *(std::max_element(gids.begin(), gids.end()));
00985   dbgOut.tprintf(1, "Vertex gids %d-%d\n", vmin, vmax);
00986 #endif
00987 
00988   // Add elements to the range passed in
00989   faces.insert(scd_box->start_element(), scd_box->start_element() + scd_box->num_elements() - 1);
00990 
00991   if (2 <= dbgOut.get_verbosity()) {
00992     assert(scd_box->boundary_complete());
00993     EntityHandle dum_ent = scd_box->start_element();
00994     rval = mbImpl->list_entities(&dum_ent, 1);
00995     ERRORR(rval, "Trouble listing first hex.");
00996 
00997     std::vector<EntityHandle> connect;
00998     rval = mbImpl->get_connectivity(&dum_ent, 1, connect);
00999     ERRORR(rval, "Trouble getting connectivity.");
01000 
01001     rval = mbImpl->list_entities(&connect[0], connect.size());
01002     ERRORR(rval, "Trouble listing element connectivity.");
01003   }
01004 
01005   Range edges;
01006   mbImpl->get_adjacencies(faces, 1, true, edges, Interface::UNION);
01007 
01008   return MB_SUCCESS;
01009 }
01010 
01011 ErrorCode ScdNCHelper::read_variables(std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
01012 {
01013   std::vector<ReadNC::VarData> vdatas;
01014   std::vector<ReadNC::VarData> vsetdatas;
01015 
01016   ErrorCode rval = read_variable_setup(var_names, tstep_nums, vdatas, vsetdatas);
01017   ERRORR(rval, "Trouble setting up read variable.");
01018 
01019   // Create COORDS tag for quads
01020   rval = create_quad_coordinate_tag();
01021   ERRORR(rval, "Trouble creating coordinate tags to entities quads");
01022 
01023   if (!vsetdatas.empty()) {
01024     rval = read_variable_to_set(vsetdatas, tstep_nums);
01025     ERRORR(rval, "Trouble read variables to set.");
01026   }
01027 
01028   if (!vdatas.empty()) {
01029     rval = read_scd_variable_to_nonset(vdatas, tstep_nums);
01030     ERRORR(rval, "Trouble read variables to entities verts/edges/faces.");
01031   }
01032 
01033   return MB_SUCCESS;
01034 }
01035 
01036 ErrorCode ScdNCHelper::read_scd_variable_to_nonset_allocate(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums)
01037 {
01038   Interface*& mbImpl = _readNC->mbImpl;
01039   std::vector<int>& dimLens = _readNC->dimLens;
01040   DebugOutput& dbgOut = _readNC->dbgOut;
01041 
01042   ErrorCode rval = MB_SUCCESS;
01043 
01044   Range* range = NULL;
01045 
01046   // Get vertices in set
01047   Range verts;
01048   rval = mbImpl->get_entities_by_dimension(_fileSet, 0, verts);
01049   ERRORR(rval, "Trouble getting vertices in set.");
01050   assert("Should only have a single vertex subrange, since they were read in one shot" &&
01051       verts.psize() == 1);
01052 
01053   Range edges;
01054   rval = mbImpl->get_entities_by_dimension(_fileSet, 1, edges);
01055   ERRORR(rval, "Trouble getting edges in set.");
01056 
01057   // Get faces in set
01058   Range faces;
01059   rval = mbImpl->get_entities_by_dimension(_fileSet, 2, faces);
01060   ERRORR(rval, "Trouble getting faces in set.");
01061   assert("Should only have a single face subrange, since they were read in one shot" &&
01062       faces.psize() == 1);
01063 
01064 #ifdef USE_MPI
01065   moab::Range faces_owned;
01066   bool& isParallel = _readNC->isParallel;
01067   if (isParallel) {
01068     ParallelComm*& myPcomm = _readNC->myPcomm;
01069     rval = myPcomm->filter_pstatus(faces, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &faces_owned);
01070     ERRORR(rval, "Trouble getting owned faces in set.");
01071   }
01072   else
01073     faces_owned = faces; // Not running in parallel, but still with MPI
01074 #endif
01075 
01076   for (unsigned int i = 0; i < vdatas.size(); i++) {
01077     // Support non-set variables with 4 dimensions like (time, lev, lat, lon)
01078     assert(4 == vdatas[i].varDims.size());
01079 
01080     // For a non-set variable, time should be the first dimension
01081     assert(tDim == vdatas[i].varDims[0]);
01082 
01083     // Set up readStarts and readCounts
01084     vdatas[i].readStarts.resize(4);
01085     vdatas[i].readCounts.resize(4);
01086 
01087     // First: time
01088     vdatas[i].readStarts[0] = 0; // This value is timestep dependent, will be set later
01089     vdatas[i].readCounts[0] = 1;
01090 
01091     // Next: lev
01092     vdatas[i].readStarts[1] = 0;
01093     vdatas[i].readCounts[1] = vdatas[i].numLev;
01094 
01095     // Finally: lat (or slat) and lon (or slon)
01096     switch (vdatas[i].entLoc) {
01097       case ReadNC::ENTLOCVERT:
01098         // Vertices
01099         vdatas[i].readStarts[2] = lDims[1];
01100         vdatas[i].readCounts[2] = lDims[4] - lDims[1] + 1;
01101         vdatas[i].readStarts[3] = lDims[0];
01102         vdatas[i].readCounts[3] = lDims[3] - lDims[0] + 1;
01103         range = &verts;
01104         break;
01105       case ReadNC::ENTLOCNSEDGE:
01106         ERRORR(MB_FAILURE, "Reading edge data not implemented yet.");
01107         break;
01108       case ReadNC::ENTLOCEWEDGE:
01109         ERRORR(MB_FAILURE, "Reading edge data not implemented yet.");
01110         break;
01111       case ReadNC::ENTLOCFACE:
01112         // Faces
01113         vdatas[i].readStarts[2] = lCDims[1];
01114         vdatas[i].readCounts[2] = lCDims[4] - lCDims[1] + 1;
01115         vdatas[i].readStarts[3] = lCDims[0];
01116         vdatas[i].readCounts[3] = lCDims[3] - lCDims[0] + 1;
01117 #ifdef USE_MPI
01118         range = &faces_owned;
01119 #else
01120         range = &faces;
01121 #endif
01122         break;
01123       case ReadNC::ENTLOCSET:
01124         // Set
01125         break;
01126       default:
01127         ERRORR(MB_FAILURE, "Unrecognized entity location type.");
01128         break;
01129     }
01130 
01131     for (unsigned int t = 0; t < tstep_nums.size(); t++) {
01132       dbgOut.tprintf(2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
01133 
01134       if (tstep_nums[t] >= dimLens[tDim]) {
01135         ERRORR(MB_INDEX_OUT_OF_RANGE, "Wrong value for a timestep number.");
01136       }
01137 
01138       // Get the tag to read into
01139       if (!vdatas[i].varTags[t]) {
01140         rval = get_tag_to_nonset(vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev);
01141         ERRORR(rval, "Trouble getting tag.");
01142       }
01143 
01144       // Get ptr to tag space
01145       void* data;
01146       int count;
01147       rval = mbImpl->tag_iterate(vdatas[i].varTags[t], range->begin(), range->end(), count, data);
01148       ERRORR(rval, "Failed to get tag iterator.");
01149       assert((unsigned)count == range->size());
01150       vdatas[i].varDatas[t] = data;
01151     }
01152 
01153     // Get variable size
01154     vdatas[i].sz = 1;
01155     for (std::size_t idx = 0; idx != vdatas[i].readCounts.size(); idx++)
01156       vdatas[i].sz *= vdatas[i].readCounts[idx];
01157   }
01158 
01159   return rval;
01160 }
01161 
01162 ErrorCode ScdNCHelper::read_scd_variable_to_nonset(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums)
01163 {
01164   DebugOutput& dbgOut = _readNC->dbgOut;
01165 
01166   ErrorCode rval = read_scd_variable_to_nonset_allocate(vdatas, tstep_nums);
01167   ERRORR(rval, "Trouble allocating read variables.");
01168 
01169   // Finally, read into that space
01170   int success;
01171   for (unsigned int i = 0; i < vdatas.size(); i++) {
01172     std::size_t sz = vdatas[i].sz;
01173 
01174     // A typical supported variable: float T(time, lev, lat, lon)
01175     // For tag values, need transpose (lev, lat, lon) to (lat, lon, lev)
01176     size_t ni = vdatas[i].readCounts[3]; // lon or slon
01177     size_t nj = vdatas[i].readCounts[2]; // lat or slat
01178     size_t nk = vdatas[i].readCounts[1]; // lev
01179 
01180     for (unsigned int t = 0; t < tstep_nums.size(); t++) {
01181       // Tag data for this timestep
01182       void* data = vdatas[i].varDatas[t];
01183 
01184       // Set readStart for each timestep along time dimension
01185       vdatas[i].readStarts[0] = tstep_nums[t];
01186 
01187       switch (vdatas[i].varDataType) {
01188         case NC_BYTE:
01189         case NC_CHAR: {
01190           std::vector<char> tmpchardata(sz);
01191           success = NCFUNCAG(_vara_text)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0], &vdatas[i].readCounts[0],
01192                                         &tmpchardata[0]);
01193           if (vdatas[i].numLev != 1)
01194             // Transpose (lev, lat, lon) to (lat, lon, lev)
01195             success = kji_to_jik(ni, nj, nk, data, &tmpchardata[0]);
01196           else {
01197             for (std::size_t idx = 0; idx != tmpchardata.size(); idx++)
01198               ((char*) data)[idx] = tmpchardata[idx];
01199           }
01200           ERRORS(success, "Failed to read char data.");
01201           break;
01202         }
01203         case NC_DOUBLE: {
01204           std::vector<double> tmpdoubledata(sz);
01205           success = NCFUNCAG(_vara_double)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0], &vdatas[i].readCounts[0],
01206                                           &tmpdoubledata[0]);
01207           if (vdatas[i].numLev != 1)
01208             // Transpose (lev, lat, lon) to (lat, lon, lev)
01209             success = kji_to_jik(ni, nj, nk, data, &tmpdoubledata[0]);
01210           else {
01211             for (std::size_t idx = 0; idx != tmpdoubledata.size(); idx++)
01212               ((double*) data)[idx] = tmpdoubledata[idx];
01213           }
01214           ERRORS(success, "Failed to read double data.");
01215           break;
01216         }
01217         case NC_FLOAT: {
01218           std::vector<float> tmpfloatdata(sz);
01219           success = NCFUNCAG(_vara_float)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0], &vdatas[i].readCounts[0],
01220                                           &tmpfloatdata[0]);
01221           if (vdatas[i].numLev != 1)
01222             // Transpose (lev, lat, lon) to (lat, lon, lev)
01223             success = kji_to_jik(ni, nj, nk, data, &tmpfloatdata[0]);
01224           else {
01225             for (std::size_t idx = 0; idx != tmpfloatdata.size(); idx++)
01226               ((float*) data)[idx] = tmpfloatdata[idx];
01227           }
01228           ERRORS(success, "Failed to read float data.");
01229           break;
01230         }
01231         case NC_INT: {
01232           std::vector<int> tmpintdata(sz);
01233           success = NCFUNCAG(_vara_int)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0], &vdatas[i].readCounts[0],
01234                                         &tmpintdata[0]);
01235           if (vdatas[i].numLev != 1)
01236             // Transpose (lev, lat, lon) to (lat, lon, lev)
01237             success = kji_to_jik(ni, nj, nk, data, &tmpintdata[0]);
01238           else {
01239             for (std::size_t idx = 0; idx != tmpintdata.size(); idx++)
01240               ((int*) data)[idx] = tmpintdata[idx];
01241           }
01242           ERRORS(success, "Failed to read int data.");
01243           break;
01244         }
01245         case NC_SHORT: {
01246           std::vector<short> tmpshortdata(sz);
01247           success = NCFUNCAG(_vara_short)(_fileId, vdatas[i].varId, &vdatas[i].readStarts[0], &vdatas[i].readCounts[0],
01248                                           &tmpshortdata[0]);
01249           if (vdatas[i].numLev != 1)
01250             // Transpose (lev, lat, lon) to (lat, lon, lev)
01251             success = kji_to_jik(ni, nj, nk, data, &tmpshortdata[0]);
01252           else {
01253             for (std::size_t idx = 0; idx != tmpshortdata.size(); idx++)
01254               ((short*) data)[idx] = tmpshortdata[idx];
01255           }
01256           ERRORS(success, "Failed to read short data.");
01257           break;
01258         }
01259         default:
01260           success = 1;
01261       }
01262 
01263       if (success)
01264         ERRORR(MB_FAILURE, "Trouble reading variable.");
01265     }
01266   }
01267 
01268   for (unsigned int i = 0; i < vdatas.size(); i++) {
01269     for (unsigned int t = 0; t < tstep_nums.size(); t++) {
01270       dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
01271       ErrorCode tmp_rval = convert_variable(vdatas[i], t);
01272       if (MB_SUCCESS != tmp_rval)
01273         rval = tmp_rval;
01274     }
01275   }
01276 
01277   // Debug output, if requested
01278   if (1 == dbgOut.get_verbosity()) {
01279     dbgOut.printf(1, "Read variables: %s", vdatas.begin()->varName.c_str());
01280     for (unsigned int i = 1; i < vdatas.size(); i++)
01281       dbgOut.printf(1, ", %s ", vdatas[i].varName.c_str());
01282     dbgOut.tprintf(1, "\n");
01283   }
01284 
01285   return rval;
01286 }
01287 
01288 ErrorCode ScdNCHelper::create_quad_coordinate_tag() {
01289   Interface*& mbImpl = _readNC->mbImpl;
01290 
01291   Range ents;
01292   ErrorCode rval = mbImpl->get_entities_by_type(_fileSet, moab::MBQUAD, ents);
01293   ERRORR(rval, "Trouble getting QUAD entity.");
01294 
01295   std::size_t numOwnedEnts = 0;
01296 #ifdef USE_MPI
01297   Range ents_owned;
01298   bool& isParallel = _readNC->isParallel;
01299   if (isParallel) {
01300     ParallelComm*& myPcomm = _readNC->myPcomm;
01301     rval = myPcomm->filter_pstatus(ents, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &ents_owned);
01302     ERRORR(rval, "Trouble getting owned QUAD entity.");
01303     numOwnedEnts = ents_owned.size();
01304   }
01305   else
01306   {
01307     numOwnedEnts = ents.size();
01308     ents_owned = ents;
01309   }
01310 #else
01311   numOwnedEnts = ents.size();
01312 #endif
01313 
01314   if (numOwnedEnts == 0)
01315     return MB_SUCCESS;
01316 
01317   assert(numOwnedEnts == ilCVals.size() * jlCVals.size());
01318   std::vector<double> coords(numOwnedEnts * 3);
01319   std::size_t pos = 0;
01320   for (std::size_t j = 0; j != jlCVals.size(); ++j) {
01321     for (std::size_t i = 0; i != ilCVals.size(); ++i) {
01322       pos = j * ilCVals.size() * 3 + i * 3;
01323       coords[pos] = ilCVals[i];
01324       coords[pos + 1] = jlCVals[j];
01325       coords[pos + 2] = 0.0;
01326     }
01327   }
01328   std::string tag_name = "COORDS";
01329   Tag tagh = 0;
01330   rval = mbImpl->tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT);
01331   ERRORR(rval, "Trouble creating COORDS tag.");
01332 
01333   void *data;
01334   int count;
01335 #ifdef USE_MPI
01336   rval = mbImpl->tag_iterate(tagh, ents_owned.begin(), ents_owned.end(), count, data);
01337 #else
01338   rval = mbImpl->tag_iterate(tagh, ents.begin(), ents.end(), count, data);
01339 #endif
01340   ERRORR(rval, "Failed to get COORDS tag iterator.");
01341   assert(count == (int)numOwnedEnts);
01342   double* quad_data = (double*) data;
01343   std::copy(coords.begin(), coords.end(), quad_data);
01344 
01345   return MB_SUCCESS;
01346 }
01347 
01348 ErrorCode UcdNCHelper::read_variables(std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
01349 {
01350   std::vector<ReadNC::VarData> vdatas;
01351   std::vector<ReadNC::VarData> vsetdatas;
01352 
01353   ErrorCode rval = read_variable_setup(var_names, tstep_nums, vdatas, vsetdatas);
01354   ERRORR(rval, "Trouble setting up read variable.");
01355 
01356   if (!vsetdatas.empty()) {
01357     rval = read_variable_to_set(vsetdatas, tstep_nums);
01358     ERRORR(rval, "Trouble read variables to set.");
01359   }
01360 
01361   if (!vdatas.empty()) {
01362 #ifdef PNETCDF_FILE
01363     // With pnetcdf support, we will use async read
01364     rval = read_ucd_variable_to_nonset_async(vdatas, tstep_nums);
01365 #else
01366     // Without pnetcdf support, we will use old read
01367     rval = read_ucd_variable_to_nonset(vdatas, tstep_nums);
01368 #endif
01369     ERRORR(rval, "Trouble read variables to entities verts/edges/faces.");
01370   }
01371 
01372   return MB_SUCCESS;
01373 }
01374 
01375 } // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines