moab
|
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