[MOAB-dev] commit/MOAB: 10 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Jun 5 11:49:41 CDT 2013
10 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/808d4e0ba385/
Changeset: 808d4e0ba385
Branch: None
User: danwu
Date: 2013-05-14 17:54:31
Summary: Refactor ReadNC class with a helper class NCHelper
Affected #: 2 files
diff --git a/src/io/ReadNC.cpp b/src/io/ReadNC.cpp
index 774f7a6..6aa41fa 100644
--- a/src/io/ReadNC.cpp
+++ b/src/io/ReadNC.cpp
@@ -41,7 +41,7 @@ ReadNC::ReadNC(Interface* impl) :
#ifdef USE_MPI
myPcomm(NULL),
#endif
- noMesh(false), noVars(false), spectralMesh(false)
+ noMesh(false), noVars(false), spectralMesh(false), helper(NULL)
{
assert(impl != NULL);
@@ -175,26 +175,17 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
if (!scdi)
return MB_FAILURE;
- // get the type of CAM file
- rval = get_nc_type(opts);
-
- if (CAM_FV == camType) {
- rval = init_FVCDscd_vals(opts, tmp_set);
- ERRORR(rval, "Trouble initializing FV grid.");
- }
- else if (CAM_SE == camType) {
- rval = init_HOMMEucd_vals();
- ERRORR(rval, "Failed to read HOMME data.");
- }
- else if (CAM_EUL == camType) {
- rval = init_EulSpcscd_vals(opts, tmp_set);
- ERRORR(rval, "Failure reading Euler grid.");
- }
- else {
- //will fill this in later for POP, CICE and CLM
- ERRORR(MB_FAILURE, "Unknown grid");
+ if (helper != NULL)
+ {
+ delete helper;
+ helper = NULL;
}
+ helper = NCHelper::get_nc_helper(fileId, this, opts);
+ rval = helper->init_nc_vals(opts, tmp_set);
+ if (MB_SUCCESS != rval)
+ return rval;
+
// Create mesh vertex/quads sequences
Range quads;
if (noMesh && !noVars) {
@@ -359,7 +350,10 @@ ErrorCode ReadNC::get_nc_type(const FileOptions &opts)
else if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
dimNames.end(), std::string("lat")) != dimNames.end()))
camType = CAM_EUL;
- else ERRORR(MB_FAILURE, "Unknown CAM grid");
+ else {
+ camType = CAM_UNKNOWN;
+ ERRORR(MB_FAILURE, "Unknown CAM grid");
+ }
}
return MB_SUCCESS;
@@ -3401,4 +3395,63 @@ ErrorCode ReadNC::create_quad_coordinate_tag(EntityHandle file_set) {
return MB_SUCCESS;
}
+NCHelper* NCHelper::get_nc_helper(int fileId, ReadNC* readNC, const FileOptions& opts)
+{
+ readNC->camType = ReadNC::NOT_CAM;
+ readNC->get_nc_type(opts);
+
+ if (ReadNC::CAM_EUL == readNC->camType)
+ return new NCHEuler(fileId, readNC);
+ else if (ReadNC::CAM_FV == readNC->camType)
+ return new NCHFV(fileId, readNC);
+ else if (ReadNC::CAM_SE == readNC->camType)
+ return new NCHHomme(fileId, readNC);
+ // will fill this in later for POP, CICE and CLM
+ else if (ReadNC::CAM_UNKNOWN == readNC->camType)
+ return new NCHUnknown(fileId, readNC);
+
+ return new NCHNotCam(fileId, readNC);
+}
+
+ErrorCode NCHEuler::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ ErrorCode rval = _readNC->init_EulSpcscd_vals(opts, file_set);
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing Euler grid.");
+
+ return rval;
+}
+
+ErrorCode NCHFV::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ ErrorCode rval = _readNC->init_FVCDscd_vals(opts, file_set);
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing FV grid.");
+
+ return rval;
+}
+
+ErrorCode NCHHomme::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ ErrorCode rval = _readNC->init_HOMMEucd_vals();
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing Homme grid.");
+
+ return rval;
+}
+
+ErrorCode NCHUnknown::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing unknown cam grid.");
+
+ return MB_FAILURE;
+}
+
+ErrorCode NCHNotCam::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing non-cam grid.");
+
+ return MB_FAILURE;
+}
+
} // namespace moab
diff --git a/src/io/ReadNC.hpp b/src/io/ReadNC.hpp
index 7bc07d5..1160315 100644
--- a/src/io/ReadNC.hpp
+++ b/src/io/ReadNC.hpp
@@ -62,12 +62,19 @@ namespace moab {
class ReadUtilIface;
class ScdInterface;
+class NCHelper;
//! Output Exodus File for VERDE
class ReadNC : public ReaderIface
{
-
+ friend class NCHelper;
+ friend class NCHEuler;
+ friend class NCHFV;
+ friend class NCHHomme;
+ friend class NCHUnknown;
+ friend class NCHNotCam;
+
public:
static ReaderIface* factory( Interface* );
@@ -416,6 +423,8 @@ private:
// read option
std::string partitionTagName;
+
+ NCHelper* helper;
};
// inline functions
@@ -424,6 +433,67 @@ inline unsigned int ReadNC::number_dimensions()
return dimVals.size();
}
+class NCHelper
+{
+public:
+ NCHelper(int fileId, ReadNC* readNC) : _fileId(fileId), _readNC(readNC) {}
+
+ static NCHelper* get_nc_helper(int fileId, ReadNC* readNC, const FileOptions& opts);
+
+ virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set) = 0;
+
+ bool is_same_file(int fileId) { return _fileId == fileId; }
+
+protected:
+ int _fileId;
+ ReadNC* _readNC;
+};
+
+class NCHEuler : public NCHelper
+{
+public:
+ NCHEuler(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+
+private:
+ virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+};
+
+class NCHFV : public NCHelper
+{
+public:
+ NCHFV(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+
+private:
+ virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+};
+
+class NCHHomme : public NCHelper
+{
+public:
+ NCHHomme(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+
+private:
+ virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+};
+
+class NCHUnknown : public NCHelper
+{
+public:
+ NCHUnknown(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+
+private:
+ virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+};
+
+class NCHNotCam : public NCHelper
+{
+public:
+ NCHNotCam(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+
+private:
+ virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+};
+
} // namespace moab
#endif
https://bitbucket.org/fathomteam/moab/commits/8ac11c91ecdd/
Changeset: 8ac11c91ecdd
Branch: None
User: danwu
Date: 2013-05-16 17:14:22
Summary: Remove get_nc_type() and embed the knowledge of how to recognize a given file type in the helper class
Affected #: 2 files
diff --git a/src/io/ReadNC.cpp b/src/io/ReadNC.cpp
index 6aa41fa..b194314 100644
--- a/src/io/ReadNC.cpp
+++ b/src/io/ReadNC.cpp
@@ -36,7 +36,7 @@ ReaderIface* ReadNC::factory(Interface* iface) {
ReadNC::ReadNC(Interface* impl) :
mbImpl(impl), CPU_WORD_SIZE(-1), IO_WORD_SIZE(-1), fileId(-1), tMin(-1), tMax(-1), iDim(-1), jDim(-1), tDim(-1), iCDim(-1),
jCDim(-1), numUnLim(-1), mCurrentMeshHandle(0), startVertex(0), startElem(0), mGlobalIdTag(0), mpFileIdTag(NULL), max_line_length(-1),
- max_str_length(-1), vertexOffset(0), dbgOut(stderr), isParallel(false), partMethod(-1), camType(NOT_CAM), isCf(false),
+ max_str_length(-1), vertexOffset(0), dbgOut(stderr), isParallel(false), partMethod(-1), isCam(false), isCf(false),
spectralOrder(-1),
#ifdef USE_MPI
myPcomm(NULL),
@@ -175,14 +175,22 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
if (!scdi)
return MB_FAILURE;
+ // Check if CF convention is being followed
+ rval = check_conventions_attribute();
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // Check if it is a CAM file
+ rval = check_source_attribute();
+ if (MB_SUCCESS != rval)
+ return rval;
+
if (helper != NULL)
- {
delete helper;
- helper = NULL;
- }
- helper = NCHelper::get_nc_helper(fileId, this, opts);
- rval = helper->init_nc_vals(opts, tmp_set);
+ helper = NCHelper::get_nc_helper(this, fileId, opts);
+
+ rval = helper->init_vals(opts, tmp_set);
if (MB_SUCCESS != rval)
return rval;
@@ -193,13 +201,13 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
ERRORR(rval, "Mesh characteristics didn't match from last read.\n");
}
else if (!noMesh) {
- if (CAM_SE == camType)
+ if (CAM_SE == helper->get_cam_type())
rval = create_ucd_verts_quads(opts, tmp_set, quads);
else
rval = create_verts_quads(scdi, tmp_set, quads);
ERRORR(rval, "Trouble creating vertices and quads.");
}
- if (noMesh && CAM_SE==camType)
+ if (noMesh && CAM_SE == helper->get_cam_type())
{
// we need to populate localGid range with the gids of vertices from the tmp_set
// localGid is important in reading the variable data into the nodes
@@ -289,76 +297,70 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
return MB_SUCCESS;
}
-
-ErrorCode ReadNC::get_nc_type(const FileOptions &opts)
+ErrorCode ReadNC::check_conventions_attribute()
{
- // check if CF convention is being followed
- std::string attname;
+ isCf = false;
+
std::map<std::string, AttData>::iterator attIt = globalAtts.find("conventions");
- if (attIt == globalAtts.end()) {
+ if (attIt == globalAtts.end())
attIt = globalAtts.find("Conventions");
- attname = std::string("Conventions");
- }
- else
- attname = std::string("conventions");
- if (attIt == globalAtts.end())
+ if (attIt == globalAtts.end()) {
ERRORR(MB_FAILURE, "File does not have conventions global attribute.\n");
+ }
unsigned int sz = attIt->second.attLen;
std::string att_data;
att_data.resize(sz + 1);
att_data[sz] = '\000';
int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
- ERRORS(success, "Failed to read attribute char data.");
+ ERRORS(success, "Failed to read conventions global attribute char data.");
if (att_data.find("CF") == std::string::npos) {
ERRORR(MB_FAILURE, "File not following known conventions.\n");
}
- else isCf = true;
-
- attIt = globalAtts.find("source");
- bool is_cam = false;
- if (attIt != globalAtts.end()) {
- sz = attIt->second.attLen;
- char* tmp_str = (char *) malloc(sz + 1);
- tmp_str[sz] = '\000';
- success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), tmp_str);
- ERRORS(success, "Failed to read attribute char data.");
- std::string tmpstr(tmp_str);
- std::string cf("CAM");
- if (tmpstr.find(cf) != std::string::npos)
- is_cam = true;
- }
- if (is_cam) {
- attIt = globalAtts.find("np");
-
- // if dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
- if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
- dimNames.end(), std::string("lat")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slon"))
- != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end()))
- camType = CAM_FV;
- // else if global attribute "np" exists then it's the HOMME grid
- else if (attIt != globalAtts.end()) {
- camType = CAM_SE;
- success = NCFUNC(get_att_int)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &spectralOrder);
- spectralOrder--; // spectral order is one less than np
- ERRORS(success, "Failed to read attribute int data.");
- if (opts.match_option("PARTITION_METHOD", "NODAL_PARTITION") == MB_SUCCESS)
- partMethod = -1;
- }
- // else if dimension names "lon" and "lat' exist then it's the Eulerian Spectral grid
- else if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
- dimNames.end(), std::string("lat")) != dimNames.end()))
- camType = CAM_EUL;
- else {
- camType = CAM_UNKNOWN;
- ERRORR(MB_FAILURE, "Unknown CAM grid");
- }
+ else
+ isCf = true;
+
+ return MB_SUCCESS;
+}
+
+ErrorCode ReadNC::check_source_attribute()
+{
+ isCam = false;
+
+ std::map<std::string, AttData>::iterator attIt = globalAtts.find("source");
+
+ if (attIt == globalAtts.end()) {
+ ERRORR(MB_FAILURE, "File does not have source global attribute.\n");
}
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ ERRORS(success, "Failed to read source global attribute char data.");
+ if (att_data.find("CAM") != std::string::npos)
+ isCam = true;
+
return MB_SUCCESS;
}
-
+
+ErrorCode ReadNC::check_np_attribute()
+{
+ std::map<std::string, AttData>::iterator attIt = globalAtts.find("np");
+ if (attIt != globalAtts.end())
+ {
+ int success = NCFUNC(get_att_int)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &spectralOrder);
+ ERRORS(success, "Failed to read np global attribute int data.");
+ spectralOrder--; // Spectral order is one less than np
+
+ return MB_SUCCESS;
+ }
+
+ return MB_FAILURE;
+}
+
ErrorCode ReadNC::load_BIL(std::string , const EntityHandle* , const FileOptions& , const Tag* ) {
/*
@@ -513,8 +515,6 @@ ErrorCode ReadNC::parse_options(const FileOptions &opts, std::vector<std::string
else
partMethod = dum;
- if (4==dum) camType = CAM_SE;
-
#endif
return MB_SUCCESS;
@@ -941,7 +941,7 @@ ErrorCode ReadNC::read_variable_setup(std::vector<std::string> &var_names, std::
std::vector<VarData> &vdatas, std::vector<VarData> &vsetdatas) {
std::map<std::string, VarData>::iterator mit;
- if (camType != CAM_SE) { // scd mesh
+ if (helper->get_cam_type() != CAM_SE) { // scd mesh
// if empty read them all
if (var_names.empty()) {
for (mit = varInfo.begin(); mit != varInfo.end(); mit++) {
@@ -1140,7 +1140,7 @@ ErrorCode ReadNC::read_variable_allocate(EntityHandle file_set, std::vector<VarD
switch (vdatas[i].entLoc) {
case 0:
// vertices
- if (camType != CAM_SE) {
+ if (helper->get_cam_type() != CAM_SE) {
// only structured mesh has j parameter that multiplies i to get total # vertices
vdatas[i].readDims[t].push_back(lDims[1]);
vdatas[i].readCounts[t].push_back(lDims[4] - lDims[1] + 1);
@@ -1249,7 +1249,7 @@ ErrorCode ReadNC::read_variables(EntityHandle file_set, std::vector<std::string>
ErrorCode rval = read_variable_setup(var_names, tstep_nums, vdatas, vsetdatas);
ERRORR(rval, "Trouble setting up read variable.");
- if (camType != CAM_SE) {
+ if (helper->get_cam_type() != CAM_SE) {
// create COORDS tag for quads
rval = create_quad_coordinate_tag(file_set);
ERRORR(rval, "Trouble creating coordinate tags to entities quads");
@@ -1262,7 +1262,7 @@ ErrorCode ReadNC::read_variables(EntityHandle file_set, std::vector<std::string>
if (!vdatas.empty()) {
#ifdef PNETCDF_FILE
- if (camType == CAM_SE) // in serial, we will use the old read, everything is contiguous
+ if (helper->get_cam_type() == CAM_SE) // in serial, we will use the old read, everything is contiguous
// in parallel, we will use async read in pnetcdf
// the other mechanism is not working, forget about it
rval = read_variable_to_nonset_async(file_set, vdatas, tstep_nums);
@@ -1462,7 +1462,7 @@ ErrorCode ReadNC::read_variable_to_nonset(EntityHandle file_set, std::vector<Var
void *data = vdatas[i].varDatas[t];
std::size_t sz = 1;
size_t ni = vdatas[i].readCounts[t][2], nj = vdatas[i].readCounts[t][3], nk = vdatas[i].readCounts[t][1];
- if (camType != CAM_SE) {
+ if (helper->get_cam_type() != CAM_SE) {
for (std::size_t idx = 0; idx != vdatas[i].readCounts[t].size(); ++idx)
sz *= vdatas[i].readCounts[t][idx];
}
@@ -1504,7 +1504,7 @@ ErrorCode ReadNC::read_variable_to_nonset(EntityHandle file_set, std::vector<Var
case NC_FLOAT: {
std::vector<float> tmpfloatdata(sz);
- if (camType != CAM_SE)
+ if (helper->get_cam_type() != CAM_SE)
{
success = NCFUNCAG(_vara_float)(fileId, vdatas[i].varId, &vdatas[i].readDims[t][0], &vdatas[i].readCounts[t][0],
&tmpfloatdata[0] NCREQ);
@@ -1795,7 +1795,7 @@ ErrorCode ReadNC::convert_variable(VarData &var_data, int tstep_num) {
void *data = var_data.varDatas[tstep_num];
std::size_t sz = 1;
- if (camType != CAM_SE) {
+ if (helper->get_cam_type() != CAM_SE) {
for (std::size_t idx = 0; idx != var_data.readCounts[tstep_num].size(); ++idx)
sz *= var_data.readCounts[tstep_num][idx];
}
@@ -3251,7 +3251,7 @@ ErrorCode ReadNC::create_tags(ScdInterface *scdi, EntityHandle file_set, const s
Tag meshTypeTag = 0;
tag_name = "__MESH_TYPE";
std::string meshTypeName;
- switch(camType)
+ switch (helper->get_cam_type())
{
case CAM_EUL: meshTypeName="CAM_EUL"; break;
case CAM_FV: meshTypeName="CAM_FV"; break;
@@ -3395,25 +3395,44 @@ ErrorCode ReadNC::create_quad_coordinate_tag(EntityHandle file_set) {
return MB_SUCCESS;
}
-NCHelper* NCHelper::get_nc_helper(int fileId, ReadNC* readNC, const FileOptions& opts)
+NCHelper* NCHelper::get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts)
{
- readNC->camType = ReadNC::NOT_CAM;
- readNC->get_nc_type(opts);
-
- if (ReadNC::CAM_EUL == readNC->camType)
- return new NCHEuler(fileId, readNC);
- else if (ReadNC::CAM_FV == readNC->camType)
- return new NCHFV(fileId, readNC);
- else if (ReadNC::CAM_SE == readNC->camType)
- return new NCHHomme(fileId, readNC);
- // will fill this in later for POP, CICE and CLM
- else if (ReadNC::CAM_UNKNOWN == readNC->camType)
- return new NCHUnknown(fileId, readNC);
-
- return new NCHNotCam(fileId, readNC);
+ // Unknown grid, will fill this in later for POP, CICE and CLM
+ if (!readNC->isCam)
+ return new NCHNotCam(readNC, fileId);
+
+ // If a CAM file, which type?
+ if (NCHEuler::can_read_file(readNC))
+ return new NCHEuler(readNC, fileId);
+ else if (NCHFV::can_read_file(readNC))
+ return new NCHFV(readNC, fileId);
+ else if (NCHHomme::can_read_file(readNC, opts))
+ return new NCHHomme(readNC, fileId);
+
+ // Unknown CAM grid
+ return new NCHUnknownCam(readNC, fileId);
}
-ErrorCode NCHEuler::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+bool NCHEuler::can_read_file(ReadNC* readNC)
+{
+ std::vector<std::string>& dimNames = readNC->dimNames;
+
+ // If dimension names "lon" AND "lat' exist then it's the Eulerian Spectral grid or the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) != dimNames.end()))
+ {
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("slon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("slat")) != dimNames.end()))
+ return false;
+ else
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHEuler::init_vals(const FileOptions& opts, EntityHandle file_set)
{
ErrorCode rval = _readNC->init_EulSpcscd_vals(opts, file_set);
if (MB_SUCCESS != rval)
@@ -3422,7 +3441,20 @@ ErrorCode NCHEuler::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
return rval;
}
-ErrorCode NCHFV::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+bool NCHFV::can_read_file(ReadNC* readNC)
+{
+ std::vector<std::string>& dimNames = readNC->dimNames;
+
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slon"))
+ != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end()))
+ return true;
+
+ return false;
+}
+
+ErrorCode NCHFV::init_vals(const FileOptions& opts, EntityHandle file_set)
{
ErrorCode rval = _readNC->init_FVCDscd_vals(opts, file_set);
if (MB_SUCCESS != rval)
@@ -3431,7 +3463,21 @@ ErrorCode NCHFV::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
return rval;
}
-ErrorCode NCHHomme::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+bool NCHHomme::can_read_file(ReadNC* readNC, const FileOptions& opts)
+{
+ // If global attribute "np" exists then it's the HOMME grid
+ ErrorCode rval = readNC->check_np_attribute();
+ if (MB_SUCCESS == rval) {
+ if (MB_SUCCESS == opts.match_option("PARTITION_METHOD", "NODAL_PARTITION"))
+ readNC->partMethod = -1;
+
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHHomme::init_vals(const FileOptions& opts, EntityHandle file_set)
{
ErrorCode rval = _readNC->init_HOMMEucd_vals();
if (MB_SUCCESS != rval)
@@ -3440,16 +3486,16 @@ ErrorCode NCHHomme::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
return rval;
}
-ErrorCode NCHUnknown::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+ErrorCode NCHUnknownCam::init_vals(const FileOptions& opts, EntityHandle file_set)
{
- _readNC->readMeshIface->report_error("%s", "Trouble initializing unknown cam grid.");
+ _readNC->readMeshIface->report_error("%s", "Unknown CAM grid.");
return MB_FAILURE;
}
-ErrorCode NCHNotCam::init_nc_vals(const FileOptions& opts, EntityHandle file_set)
+ErrorCode NCHNotCam::init_vals(const FileOptions& opts, EntityHandle file_set)
{
- _readNC->readMeshIface->report_error("%s", "Trouble initializing non-cam grid.");
+ _readNC->readMeshIface->report_error("%s", "Unknown grid.");
return MB_FAILURE;
}
diff --git a/src/io/ReadNC.hpp b/src/io/ReadNC.hpp
index 1160315..cf6fb7d 100644
--- a/src/io/ReadNC.hpp
+++ b/src/io/ReadNC.hpp
@@ -72,7 +72,7 @@ class ReadNC : public ReaderIface
friend class NCHEuler;
friend class NCHFV;
friend class NCHHomme;
- friend class NCHUnknown;
+ friend class NCHUnknownCam;
friend class NCHNotCam;
public:
@@ -259,10 +259,13 @@ private:
ErrorCode get_BIL_dir();
bool BIL_mode_enabled(const char * file_name);
-
- //! parse various options and variables and attributes to infer CAM file type
- ErrorCode get_nc_type(const FileOptions &opts);
-
+
+ ErrorCode check_conventions_attribute();
+
+ ErrorCode check_source_attribute();
+
+ ErrorCode check_np_attribute();
+
template <typename T> ErrorCode kji_to_jik(size_t ni, size_t nj, size_t nk, void *dest, T *source)
{
size_t nik = ni * nk, nij = ni * nj;
@@ -387,8 +390,8 @@ private:
//! partitioning method
int partMethod;
- //! if a CAM file, which type
- CamType camType;
+ //! true if a CAM file
+ bool isCam;
//! true if file is marked as following CF metadata conventions
bool isCf;
@@ -433,65 +436,81 @@ inline unsigned int ReadNC::number_dimensions()
return dimVals.size();
}
+// Helper class to isolate reading of several different nc file formats
class NCHelper
{
public:
- NCHelper(int fileId, ReadNC* readNC) : _fileId(fileId), _readNC(readNC) {}
+ NCHelper(ReadNC* readNC, int fileId) : _readNC(readNC), _fileId(fileId) {}
- static NCHelper* get_nc_helper(int fileId, ReadNC* readNC, const FileOptions& opts);
+ static NCHelper* get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts);
- virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set) = 0;
-
- bool is_same_file(int fileId) { return _fileId == fileId; }
+ virtual ReadNC::CamType get_cam_type() = 0;
+ virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set) = 0;
protected:
- int _fileId;
ReadNC* _readNC;
+ int _fileId;
};
+// Child helper class for Eulerian Spectral grid (CAM_EUL)
class NCHEuler : public NCHelper
{
public:
- NCHEuler(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+ NCHEuler(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+ static bool can_read_file(ReadNC* readNC);
private:
- virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+ virtual ReadNC::CamType get_cam_type() { return ReadNC::CAM_EUL; }
+ virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
};
+// Child helper class for Finite Volume grid (CAM_FV)
class NCHFV : public NCHelper
{
public:
- NCHFV(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+ NCHFV(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+ static bool can_read_file(ReadNC* readNC);
private:
- virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+ virtual ReadNC::CamType get_cam_type() { return ReadNC::CAM_FV; }
+ virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
};
+// Child helper class for HOMME grid (CAM_SE)
class NCHHomme : public NCHelper
{
public:
- NCHHomme(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+ NCHHomme(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+ static bool can_read_file(ReadNC* readNC, const FileOptions& opts);
private:
- virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+ virtual ReadNC::CamType get_cam_type() { return ReadNC::CAM_SE; }
+ virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
};
-class NCHUnknown : public NCHelper
+// Child helper class for unknown CAM grid (CAM_UNKNOWN)
+class NCHUnknownCam : public NCHelper
{
public:
- NCHUnknown(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+ NCHUnknownCam(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
private:
- virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+ virtual ReadNC::CamType get_cam_type() { return ReadNC::CAM_UNKNOWN; }
+ virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
};
+// Child helper class for unkown grid (NOT_CAM)
class NCHNotCam : public NCHelper
{
public:
- NCHNotCam(int fileId, ReadNC* readNC) : NCHelper(fileId, readNC) {}
+ NCHNotCam(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
private:
- virtual ErrorCode init_nc_vals(const FileOptions& opts, EntityHandle file_set);
+ virtual ReadNC::CamType get_cam_type() { return ReadNC::NOT_CAM; }
+ virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
};
} // namespace moab
https://bitbucket.org/fathomteam/moab/commits/c12562f1e6a2/
Changeset: c12562f1e6a2
Branch: None
User: danwu
Date: 2013-05-16 19:51:55
Summary: Merged fathomteam/moab into master
Affected #: 4 files
diff --git a/src/AEntityFactory.cpp b/src/AEntityFactory.cpp
index 3d04244..577162e 100644
--- a/src/AEntityFactory.cpp
+++ b/src/AEntityFactory.cpp
@@ -788,7 +788,7 @@ ErrorCode AEntityFactory::get_down_adjacency_elements(EntityHandle source_entity
ErrorCode AEntityFactory::get_down_adjacency_elements_poly(EntityHandle source_entity,
const unsigned int target_dimension,
std::vector<EntityHandle> &target_entities,
- const bool /*create_if_missing*/,
+ const bool create_if_missing,
const int /*create_adjacency_option*/)
{
@@ -843,6 +843,34 @@ ErrorCode AEntityFactory::get_down_adjacency_elements_poly(EntityHandle source_e
target_entities.push_back(*adj_edges.begin());
}
}
+ else
+ {
+ // we have no adjacent edge yet; we need to create one and also add
+ // them to the adjacency of the vertices
+ if (create_if_missing)
+ {
+ EntityHandle newEdge;
+ EntityHandle v[2] = {vertex_array[i], vertex_array[i+1]};
+ result = thisMB->create_element(MBEDGE, v, 2,
+ newEdge);
+ if (MB_SUCCESS!=result)
+ return result;
+ // we also need to add explicit adjacency, so next time we do not
+ // create again (because we do not find the edge if it is not adjacent to the
+ // vertices
+ // if (create_adjacency_option >= 0)
+ //{
+ result = add_adjacency(v[0],newEdge);
+ if (MB_SUCCESS!=result)
+ return result;
+ result = add_adjacency(v[1],newEdge);
+ if (MB_SUCCESS!=result)
+ return result;
+ target_entities.push_back(newEdge);
+ //}
+
+ }
+ }
}
return result;
}
@@ -855,7 +883,7 @@ ErrorCode AEntityFactory::get_down_adjacency_elements_poly(EntityHandle source_e
std::vector<EntityHandle> dum_vec;
result = thisMB->get_connectivity(&source_entity, 1, dum_vec);
if (MB_SUCCESS != result) return result;
- result = thisMB->get_adjacencies(&dum_vec[0], dum_vec.size(), 1, false,
+ result = thisMB->get_adjacencies(&dum_vec[0], dum_vec.size(), 1, create_if_missing,
target_entities, Interface::UNION);
return result;
}
diff --git a/test/MBTest.cpp b/test/MBTest.cpp
index e695f99..7394c63 100644
--- a/test/MBTest.cpp
+++ b/test/MBTest.cpp
@@ -5950,6 +5950,97 @@ ErrorCode mb_poly_adjacency_test()
return moab.check_adjacencies();
}
+ErrorCode mb_poly_adjacency_test2()
+{
+ ErrorCode rval;
+ Core moab;
+ Interface *mbImpl = &moab;
+
+ // make a polyhedra in shape of a cube with a corner cut
+
+ /*
+ *
+ * 7 ------------- 8
+ * . | . |
+ * . . | initial cube had 8 vertices; box [000 - 222]
+ * 5 ------------- 6 | cut a triangle 123 on lower corner
+ * | | | |
+ * | | | faces of the polyhedra (pointing outward)
+ * | | | | (123) (24653) (4 10 8 6)
+ * 3 | | (9 10 4 2 1) (7 8 10 9) (5687)
+ * \ | | | (1 3 5 7 9)
+ * 9 - - -| - 10
+ * \ . | .
+ * 1, | .
+ * (000) '2 ---- 4
+ *
+ */
+ double coords[] = {0, 1, 0, // vertex 1
+ 1, 0, 0, // vertex 2
+ 0, 0 ,1, // vertex 3
+ 2, 0, 0,
+ 0, 0, 2,
+ 2, 0, 2, // vertex 6
+ 0, 2, 2,
+ 2, 2, 2, // vertex 8
+ 0, 2, 0, // vertex 9
+ 2, 2, 0 // vertex 9
+
+ };
+ EntityHandle verts[10], polygons[7], polyhedron;
+
+ for (int i = 0; i < 10; i++) {
+ rval = mbImpl->create_vertex(&coords[3*i], verts[i]);
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+
+ EntityHandle connect[]= {1, 2, 3, // poly 1
+ 2, 4, 6, 5, 3, // poly 2
+ 4, 10, 8, 6, // polygon 3...
+ 9, 10, 4, 2, 1,
+ 7, 8, 10, 9,
+ 5, 6, 8, 7,
+ 1, 3, 5, 7, 9 // polygon 7
+ }; // we know the handles directly
+
+ int num_verts[7]={3, 5, 4, 5, 4, 4, 5};
+ int start_indx[7];
+ start_indx[0]= 0;
+ for (int i=0; i<6; i++)
+ start_indx[i+1]=start_indx[i]+num_verts[i];
+ for (int j = 0; j < 7; j++) {
+ rval = mbImpl->create_element(MBPOLYGON, &connect[start_indx[j]],
+ num_verts[j], polygons[j]);
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+ rval = mbImpl->create_element(MBPOLYHEDRON, polygons, 7, polyhedron);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // create the aentities
+ Range dum_range;
+ for (int dim = 0; dim < 3; dim++) {
+ dum_range.clear();
+ rval = mbImpl->get_adjacencies(&polyhedron, 1, dim, true, dum_range, Interface::UNION);
+ if (MB_SUCCESS != rval)
+ return rval;
+ std::cout << "\n dimension:" << dim << " " << dum_range.size() << " entities";
+ }
+ std::cout << "\n";
+ /*rval=mbImpl->write_mesh("polyhedra.vtk");
+ if (MB_SUCCESS != rval)
+ return rval;*/
+ // delete the polyhedron
+ rval = mbImpl->delete_entities(&polyhedron, 1);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // check adjacencies
+ return moab.check_adjacencies();
+}
+
ErrorCode mb_memory_use_test()
{
Core mb;
@@ -8152,6 +8243,7 @@ int main(int argc, char* argv[])
RUN_TEST( mb_split_test );
RUN_TEST( mb_range_seq_intersect_test );
RUN_TEST( mb_poly_adjacency_test );
+ RUN_TEST( mb_poly_adjacency_test2 );
RUN_TEST( mb_memory_use_test );
RUN_TEST( mb_skin_curve_test );
RUN_TEST( mb_skin_curve_adj_test );
diff --git a/tools/mbcslam/Makefile.am b/tools/mbcslam/Makefile.am
index 2d37b4a..f7785e7 100644
--- a/tools/mbcslam/Makefile.am
+++ b/tools/mbcslam/Makefile.am
@@ -40,7 +40,10 @@ cfgdir = $(libdir)
TESTS = intx_on_sphere_test intx_in_plane_test spec_visu_test spherical_area_test \
case1_test
-noinst_PROGRAMS = process_arm cslam_par_test
+noinst_PROGRAMS = cslam_par_test
+if NETCDF_FILE
+ noinst_PROGRAMS += process_arm
+endif
check_PROGRAMS = $(TESTS)
intx_on_sphere_test_SOURCES = intx_on_sphere_test.cpp
diff --git a/tools/mbcslam/intx_in_plane_test.cpp b/tools/mbcslam/intx_in_plane_test.cpp
index 774e583..d91ea79 100644
--- a/tools/mbcslam/intx_in_plane_test.cpp
+++ b/tools/mbcslam/intx_in_plane_test.cpp
@@ -25,6 +25,7 @@ int main(int argc, char* argv[])
const char *filename_mesh1 = STRINGIFY(SRCDIR) "/m1.vtk";
const char *filename_mesh2 = STRINGIFY(SRCDIR) "/m2.vtk";
const char *newFile = "intx1.vtk";
+ const char *edgesFile = "polyWithEdges.vtk";
if (argc == 4)
{
filename_mesh1 = argv[1];
@@ -86,6 +87,13 @@ int main(int argc, char* argv[])
return 1;
std::cout << "number of edges:" << edges.size() << "\n";
+ // add edges to the output set
+ rval = mb->add_entities(outputSet, edges);
+ if (MB_SUCCESS != rval)
+ return 1;
+ rval = mb->write_mesh(edgesFile, &outputSet, 1);
+ if (MB_SUCCESS != rval)
+ return 1;
return 0;
https://bitbucket.org/fathomteam/moab/commits/610a825c4f22/
Changeset: 610a825c4f22
Branch: None
User: danwu
Date: 2013-05-20 16:23:38
Summary: Merged fathomteam/moab into master
Affected #: 3 files
diff --git a/itaps/imesh/FindAdjacencyF90.F90 b/itaps/imesh/FindAdjacencyF90.F90
index 78578a7..8115ef1 100644
--- a/itaps/imesh/FindAdjacencyF90.F90
+++ b/itaps/imesh/FindAdjacencyF90.F90
@@ -69,7 +69,7 @@ program findadjacency
vert_uses = vert_uses + iverts_size
- if (iverts_size .ne. 0) call free(rpverts)
+ if (iverts_size .ne. 0) call iMesh_freeMemory(%VAL(mesh), rpverts)
end do
! now get adjacencies in one big block
@@ -81,9 +81,9 @@ program findadjacency
offsets_size, ierr)
CHECK("Failure in getEntArrAdj")
- if (allverts_size .ne. 0) call free(rpallverts);
- if (offsets_size .ne. 0) call free(ipoffsets);
- if (ents_size .ne. 0) call free(rpents);
+ if (allverts_size .ne. 0) call iMesh_freeMemory(%VAL(mesh), rpallverts);
+ if (offsets_size .ne. 0) call iMesh_freeMemory(%VAL(mesh), ipoffsets);
+ if (ents_size .ne. 0) call iMesh_freeMemory(%VAL(mesh), rpents);
! compare results of two calling methods
if (allverts_size .ne. vert_uses) then
@@ -119,6 +119,9 @@ program findadjacency
!
write(*, *) "num coords: ", ents_size, " few coords: ", (coords(i), i=0, ents_size/100)
+ call iMesh_freeMemory(%VAL(mesh), verths);
+ call iMesh_freeMemory(%VAL(mesh), pcoord);
+
call iMesh_dtor(%VAL(mesh), ierr)
CHECK("Failed to destroy interface")
diff --git a/itaps/imesh/iMesh_MOAB.cpp b/itaps/imesh/iMesh_MOAB.cpp
index 963e740..6d55492 100644
--- a/itaps/imesh/iMesh_MOAB.cpp
+++ b/itaps/imesh/iMesh_MOAB.cpp
@@ -3423,3 +3423,13 @@ void iMesh_createStructuredMesh(iMesh_Instance instance,
RETURN(iBase_SUCCESS);
}
+void iMesh_freeMemory(
+ iMesh_Instance instance,
+ /**< [in] iMesh instance handle */
+ void ** ptrToMem)
+{
+ free(*ptrToMem);
+ *ptrToMem=0;
+ return;
+}
+
diff --git a/itaps/imesh/iMesh_extensions.h b/itaps/imesh/iMesh_extensions.h
index 8be6d22..c85571c 100644
--- a/itaps/imesh/iMesh_extensions.h
+++ b/itaps/imesh/iMesh_extensions.h
@@ -399,6 +399,15 @@ void iMesh_createStructuredMesh(
int *err
/**< [out] Error flag. */
);
+/***************************************************************************//**
+ * \brief Free memory allocated with malloc
+ *
+ ******************************************************************************/
+
+void iMesh_freeMemory(
+ iMesh_Instance instance,
+ /**< [in] iMesh instance handle */
+ void ** ptrToMem);
/***************************************************************************//**
* \defgroup ScdMesh Structured Mesh
https://bitbucket.org/fathomteam/moab/commits/2b13269704c3/
Changeset: 2b13269704c3
Branch: None
User: danwu
Date: 2013-05-24 15:42:24
Summary: Merged fathomteam/moab into master
Affected #: 5 files
diff --git a/config/compiler.m4 b/config/compiler.m4
index 18df411..81596a6 100644
--- a/config/compiler.m4
+++ b/config/compiler.m4
@@ -126,8 +126,6 @@ fi
if test "xno" != "x$CHECK_FC"; then
AC_PROG_FC
AC_PROG_F77
- AC_F77_LIBRARY_LDFLAGS
- AC_FC_LIBRARY_LDFLAGS
fi
]) # FATHOM_CHECK_COMPILERS
@@ -523,6 +521,12 @@ case "$cc_compiler:$host_cpu" in
case "$target_vendor" in
bgp)
FATHOM_CC_32BIT=-q32
+ FATHOM_CC_64BIT=-q64
+ AR="ar"
+ NM="nm -B"
+ ;;
+ bgq)
+ FATHOM_CC_32BIT=-q32
FATHOM_CC_64BIT=-q64
AR="ar"
NM="nm -B"
diff --git a/configure.ac b/configure.ac
index 1524b06..c9d6387 100644
--- a/configure.ac
+++ b/configure.ac
@@ -67,14 +67,6 @@ if test "xyes" = "x$ENABLE_FORTRAN" && test "x" != "x$FC"; then
AC_FC_WRAPPERS
fi
-if test "xyes" = "x$ENABLE_FORTRAN" && test "x" != "x$FLIBS"; then
- LIBS="$LIBS $FLIBS"
-fi
-
-if test "xyes" = "x$ENABLE_FORTRAN" && test "x" != "x$FCLIBS"; then
- LIBS="$LIBS $FCLIBS"
-fi
-
################################################################################
# Check for need for extra flags to support cray pointers
################################################################################
@@ -1230,10 +1222,10 @@ AC_OUTPUT
AC_MSG_RESULT([C: $CC $CFLAGS $CPPFLAGS])
AC_MSG_RESULT([C++: $CXX $CXXFLAGS $CPPFLAGS])
if test "x" != "x$FC"; then
- AC_MSG_RESULT([Fortran: $FC $FCFLAGS])
+ AC_MSG_RESULT([Fortran90: $FC $FCFLAGS])
fi
if test "x" != "x$F77"; then
- AC_MSG_RESULT([Fortran: $F77 $FFLAGS])
+ AC_MSG_RESULT([Fortran77: $F77 $FFLAGS])
fi
if test "x$HAVE_HDF5" = "xno"; then
diff --git a/doc/MOAB-UG.doc b/doc/MOAB-UG.doc
index 567a315..281d619 100644
Binary files a/doc/MOAB-UG.doc and b/doc/MOAB-UG.doc differ
diff --git a/itaps/imesh/Makefile.am b/itaps/imesh/Makefile.am
index 97c76bb..647b8c1 100644
--- a/itaps/imesh/Makefile.am
+++ b/itaps/imesh/Makefile.am
@@ -45,9 +45,6 @@ if PARALLEL
# ftest_DEPENDENCIES = libiMesh.la $(top_builddir)/libMOAB.la
endif
-FCLINK = $(CXXLINK)
-F77LINK = $(CXXLINK)
-
TESTS = $(check_PROGRAMS)
LDADD = libiMesh.la $(top_builddir)/src/libMOAB.la ${MOAB_CXX_LINKFLAGS} ${MOAB_CXX_LIBS}
TESTDEPS = libiMesh.la $(top_builddir)/src/libMOAB.la
diff --git a/src/io/Tqdcfr.cpp b/src/io/Tqdcfr.cpp
index 02a4bd9..ef01b16 100644
--- a/src/io/Tqdcfr.cpp
+++ b/src/io/Tqdcfr.cpp
@@ -18,6 +18,8 @@
#include "moab/Range.hpp"
#include "FileOptions.hpp"
#include <iostream>
+#include <string>
+
#ifdef USE_MPI
#include "moab_mpi.h"
@@ -542,33 +544,58 @@ ErrorCode Tqdcfr::read_nodeset(const unsigned int nsindex,
ns_entities, excl_entities);
if (MB_SUCCESS != result) return result;
}
-
// check for more data
+ // <- likask
if (num_read < nodeseth->nsLength) {
FREADC(2); num_read += 2;
if (char_buf[0] == 'i' && char_buf[1] == 'd') {
FREADI(1); num_read += sizeof(int);
//uid = int_buf[0];
- }
-
- if (num_read < nodeseth->nsLength) {
- // check for bc_data
- FREADC(2); num_read += 2;
+ } else {
if (char_buf[0] == 'b' && char_buf[1] == 'c') {
- FREADI(1); num_read += sizeof(int);
- int num_bcs = int_buf[0];
+ FREADI(1); num_read += sizeof(int);
+ int num_bcs = uint_buf[0];
bc_data.resize(num_bcs);
FREADCA(num_bcs, &bc_data[0]); num_read += num_bcs;
}
}
}
+ if(debug) { // <-likask
+ nodeseth->print();
+ if(!bc_data.empty()) {
+ std::cout << "bc_data = ";
+ std::vector<char>::iterator vit = bc_data.begin();
+ for(;vit!=bc_data.end();vit++) {
+ std::cout << std::hex << (int)((unsigned char)*vit) << " ";
+ }
+ std::cout << ": ";
+ vit = bc_data.begin();
+ for(;vit!=bc_data.end();vit++) {
+ std::cout << *vit;
+ }
+ std::cout << std::endl;
+ }
+ } // <-likask
+
// and put entities into this nodeset's set
ErrorCode result = put_into_set(nodeseth->setHandle, ns_entities, excl_entities);
if (MB_SUCCESS != result) return result;
result = get_names(model->nodesetMD, nsindex, nodeseth->setHandle);
if (MB_SUCCESS != result) return result;
+
+ // <-likask
+ const int def_bc_data_len = 0;
+ std::string tag_name = std::string(DIRICHLET_SET_TAG_NAME)+"__BC_DATA";
+ Tag nbc_data;
+ result = mdbImpl->tag_get_handle(tag_name.c_str(),def_bc_data_len,MB_TYPE_OPAQUE,
+ nbc_data,MB_TAG_CREAT|MB_TAG_SPARSE|MB_TAG_BYTES|MB_TAG_VARLEN,NULL);
+ if (MB_SUCCESS != result) return result;
+ void const* tag_data[] = { &(bc_data[0]) };
+ int tag_size = bc_data.size();
+ result = mdbImpl->tag_set_by_ptr(nbc_data,&nodeseth->setHandle,1,tag_data,&tag_size);
+ if (MB_SUCCESS != result) return result;
return result;
}
@@ -676,22 +703,49 @@ ErrorCode Tqdcfr::read_sideset(const unsigned int ssindex,
if (char_buf[0] == 'i' && char_buf[1] == 'd') {
FREADI(1); num_read += sizeof(int);
//uid = int_buf[0];
- }
-
- if (num_read < sideseth->ssLength) {
+ } else {
// check for bc_data
- FREADC(2); num_read += 2;
if (char_buf[0] == 'b' && char_buf[1] == 'c') {
FREADI(1); num_read += sizeof(int);
- int num_bcs = int_buf[0];
+ int num_bcs = uint_buf[0];
bc_data.resize(num_bcs);
FREADCA(num_bcs, &bc_data[0]); num_read += num_bcs;
}
}
}
+ if(debug) { // <-likask
+ sideseth->print();
+ if(!bc_data.empty()) {
+ std::cout << "bc_data = ";
+ std::vector<char>::iterator vit = bc_data.begin();
+ for(;vit!=bc_data.end();vit++) {
+ std::cout << std::hex << (int)((unsigned char)*vit) << " ";
+ }
+ std::cout << ": ";
+ vit = bc_data.begin();
+ for(;vit!=bc_data.end();vit++) {
+ std::cout << *vit;
+ }
+ std::cout << std::endl;
+ }
+ } // <-likask
+
result = get_names(model->sidesetMD, ssindex, sideseth->setHandle);
if (MB_SUCCESS != result) return result;
+
+ // <-likask
+ const int def_bc_data_len = 0;
+ std::string tag_name = std::string(NEUMANN_SET_TAG_NAME)+"__BC_DATA";
+ Tag nbc_data;
+ result = mdbImpl->tag_get_handle(tag_name.c_str(),def_bc_data_len,MB_TYPE_OPAQUE,
+ nbc_data,MB_TAG_CREAT|MB_TAG_SPARSE|MB_TAG_BYTES|MB_TAG_VARLEN,NULL);
+ if (MB_SUCCESS != result) return result;
+ void const* tag_data[] = { &(bc_data[0]) };
+ int tag_size = bc_data.size();
+ result = mdbImpl->tag_set_by_ptr(nbc_data,&sideseth->setHandle,1,tag_data,&tag_size);
+ if (MB_SUCCESS != result) return result;
+
return MB_SUCCESS;
}
@@ -876,17 +930,6 @@ ErrorCode Tqdcfr::read_block(const unsigned int blindex,
FREADI(1); num_read += sizeof(int);
//uid = int_buf[0];
}
-
- if (num_read < blockh->blockLength) {
- // check for bc_data
- FREADC(2); num_read += 2;
- if (char_buf[0] == 'b' && char_buf[1] == 'c') {
- FREADI(1); num_read += sizeof(int);
- int num_bcs = int_buf[0];
- bc_data.resize(num_bcs);
- FREADCA(num_bcs, &bc_data[0]); num_read += num_bcs;
- }
- }
}
result = get_names(model->blockMD, blindex, blockh->setHandle);
@@ -1928,6 +1971,7 @@ ErrorCode Tqdcfr::NodesetHeader::read_info_header(const unsigned int model_offse
&(nodeset_headers[i].setHandle), 1,
dirichlet_category);
if (MB_SUCCESS != result) return result;
+
}
@@ -2674,7 +2718,7 @@ void Tqdcfr::BlockHeader::print()
std::cout << "blockDim = " << blockDim << std::endl;
std::cout << "setHandle = " << setHandle << std::endl;
std::cout << "blockEntityType = " << blockEntityType << std::endl;
- }
+}
Tqdcfr::NodesetHeader::NodesetHeader()
: nsID(0), memCt(0), memOffset(0), memTypeCt(0), pointSym(0), nsCol(0), nsLength(0),
@@ -2691,7 +2735,7 @@ void Tqdcfr::NodesetHeader::print()
std::cout << "nsCol = " << nsCol << std::endl;
std::cout << "nsLength = " << nsLength << std::endl;
std::cout << "setHandle = " << setHandle << std::endl;
- }
+}
Tqdcfr::SidesetHeader::SidesetHeader()
: ssID(0), memCt(0), memOffset(0), memTypeCt(0), numDF(0), ssCol(0), useShell(0), ssLength(0),
@@ -2797,7 +2841,7 @@ using namespace moab;
int main(int argc, char* argv[])
{
#ifdef USE_MPI
- int err = MPI_Init(&argc, &argv);
+ MPI_Init(&argc, &argv);
#endif
// Check command line arg
const char* file = STRINGIFY(SRCDIR) "/brick_cubit10.2.cub";
@@ -2819,9 +2863,9 @@ int main(int argc, char* argv[])
std::cout << "Success." << std::endl;
else {
std::cout << "load_file returned error:" << std::endl;
- std::string err;
- result = my_impl->get_last_error(err);
- if (MB_SUCCESS == result) std::cout << err << std::endl;
+ std::string errs;
+ result = my_impl->get_last_error(errs);
+ if (MB_SUCCESS == result) std::cout << errs << std::endl;
else std::cout << "(no message)" << std::endl;
}
https://bitbucket.org/fathomteam/moab/commits/cc9d96a72062/
Changeset: cc9d96a72062
Branch: None
User: danwu
Date: 2013-05-29 16:14:30
Summary: Move NCHelper class to separate files. Update ReadNC and NCHelper based on recent code review
Affected #: 5 files
diff --git a/src/io/Makefile.am b/src/io/Makefile.am
index 9c76524..5a0a6f3 100644
--- a/src/io/Makefile.am
+++ b/src/io/Makefile.am
@@ -16,8 +16,9 @@ if NETCDF_FILE
MOAB_NETCDF_SRCS = ReadNCDF.cpp ReadNCDF.hpp \
WriteNCDF.cpp WriteNCDF.hpp \
WriteSLAC.cpp WriteSLAC.hpp \
- ReadNC.cpp ReadNC.hpp \
- ReadGCRM.cpp ReadGCRM.hpp
+ ReadNC.cpp ReadNC.hpp \
+ NCHelper.cpp NCHelper.hpp \
+ ReadGCRM.cpp ReadGCRM.hpp
else
MOAB_NETCDF_SRCS =
endif
@@ -25,7 +26,8 @@ endif
if PNETCDF_FILE
if !NETCDF_FILE
MOAB_NETCDF_SRCS += ReadNC.cpp ReadNC.hpp \
- ReadGCRM.cpp ReadGCRM.hpp
+ NCHelper.cpp NCHelper.hpp \
+ ReadGCRM.cpp ReadGCRM.hpp
endif
endif
diff --git a/src/io/NCHelper.cpp b/src/io/NCHelper.cpp
new file mode 100644
index 0000000..d95cd34
--- /dev/null
+++ b/src/io/NCHelper.cpp
@@ -0,0 +1,112 @@
+#include "NCHelper.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+
+namespace moab {
+
+NCHelper* NCHelper::get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts)
+{
+ // Not a CAM file, will fill this in later for POP, CICE and CLM
+ if (!readNC->isCam)
+ return new (std::nothrow) NCHNotCam(readNC, fileId);
+
+ // If a CAM file, which type?
+ if (NCHEuler::can_read_file(readNC))
+ return new (std::nothrow) NCHEuler(readNC, fileId);
+ else if (NCHFV::can_read_file(readNC))
+ return new (std::nothrow) NCHFV(readNC, fileId);
+ else if (NCHHomme::can_read_file(readNC, opts))
+ return new (std::nothrow) NCHHomme(readNC, fileId);
+
+ // Unknown CAM grid
+ return new (std::nothrow) NCHUnknownCam(readNC, fileId);
+}
+
+bool NCHEuler::can_read_file(ReadNC* readNC)
+{
+ std::vector<std::string>& dimNames = readNC->dimNames;
+
+ // If dimension names "lon" AND "lat' exist then it's the Eulerian Spectral grid or the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) != dimNames.end()))
+ {
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("slon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("slat")) != dimNames.end()))
+ return false;
+ else
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHEuler::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ ErrorCode rval = _readNC->init_EulSpcscd_vals(opts, file_set);
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing Euler grid.");
+
+ return rval;
+}
+
+bool NCHFV::can_read_file(ReadNC* readNC)
+{
+ std::vector<std::string>& dimNames = readNC->dimNames;
+
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slon"))
+ != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end()))
+ return true;
+
+ return false;
+}
+
+ErrorCode NCHFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ ErrorCode rval = _readNC->init_FVCDscd_vals(opts, file_set);
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing FV grid.");
+
+ return rval;
+}
+
+bool NCHHomme::can_read_file(ReadNC* readNC, const FileOptions& opts)
+{
+ // If global attribute "np" exists then it's the HOMME grid
+ ErrorCode rval = readNC->check_np_attribute();
+ if (MB_SUCCESS == rval) {
+ if (MB_SUCCESS == opts.match_option("PARTITION_METHOD", "NODAL_PARTITION"))
+ readNC->partMethod = -1;
+
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHHomme::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ ErrorCode rval = _readNC->init_HOMMEucd_vals();
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing Homme grid.");
+
+ return rval;
+}
+
+ErrorCode NCHUnknownCam::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ _readNC->readMeshIface->report_error("%s", "Unknown CAM grid.");
+
+ return MB_FAILURE;
+}
+
+ErrorCode NCHNotCam::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ _readNC->readMeshIface->report_error("%s", "Unknown grid.");
+
+ return MB_FAILURE;
+}
+
+} // namespace moab
diff --git a/src/io/NCHelper.hpp b/src/io/NCHelper.hpp
new file mode 100644
index 0000000..10d6cdf
--- /dev/null
+++ b/src/io/NCHelper.hpp
@@ -0,0 +1,158 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelper.hpp
+//
+// Purpose : Climate NC file helper
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPER_HPP
+#define NCHELPER_HPP
+
+#include "ReadNC.hpp"
+
+namespace moab {
+
+//! Helper class to isolate reading of several different nc file formats
+class NCHelper
+{
+public:
+ NCHelper(ReadNC* readNC, int fileId) : _readNC(readNC), _fileId(fileId) {}
+
+ static NCHelper* get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts);
+
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set) = 0;
+
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads) = 0;
+
+ virtual ErrorCode init_local_gid(EntityHandle file_set) = 0;
+
+ virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums) = 0;
+
+ virtual std::string get_mesh_type_name() = 0;
+
+protected:
+ ReadNC* _readNC;
+
+ int _fileId;
+};
+
+//! Child helper class for Eulerian Spectral grid (CAM_EUL)
+class NCHEuler : public NCHelper
+{
+public:
+ NCHEuler(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+ static bool can_read_file(ReadNC* readNC);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+ { return _readNC->create_scd_verts_quads(scdi, file_set, quads); }
+
+ virtual ErrorCode init_local_gid(EntityHandle file_set)
+ { return MB_SUCCESS; }
+
+ virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
+ { return _readNC->read_variables(file_set, var_names, tstep_nums, true); }
+
+ virtual std::string get_mesh_type_name()
+ { return "CAM_EUL"; }
+};
+
+//! Child helper class for Finite Volume grid (CAM_FV)
+class NCHFV : public NCHelper
+{
+public:
+ NCHFV(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+ static bool can_read_file(ReadNC* readNC);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+ { return _readNC->create_scd_verts_quads(scdi, file_set, quads); }
+
+ virtual ErrorCode init_local_gid(EntityHandle file_set)
+ { return MB_SUCCESS; }
+
+ virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
+ { return _readNC->read_variables(file_set, var_names, tstep_nums, true); }
+
+ virtual std::string get_mesh_type_name()
+ { return "CAM_FV"; }
+};
+
+//! Child helper class for HOMME grid (CAM_SE)
+class NCHHomme : public NCHelper
+{
+public:
+ NCHHomme(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+ static bool can_read_file(ReadNC* readNC, const FileOptions& opts);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+ { return _readNC->create_ucd_verts_quads(opts, file_set, quads); }
+
+ virtual ErrorCode init_local_gid(EntityHandle file_set)
+ { return _readNC->init_local_gid(file_set); }
+
+ virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
+ { return _readNC->read_variables(file_set, var_names, tstep_nums, false); }
+
+ virtual std::string get_mesh_type_name()
+ { return "CAM_SE"; }
+};
+
+//! Child helper class for unknown CAM grid (CAM_UNKNOWN)
+class NCHUnknownCam : public NCHelper
+{
+public:
+ NCHUnknownCam(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+ { return MB_FAILURE; }
+
+ virtual ErrorCode init_local_gid(EntityHandle file_set)
+ { return MB_FAILURE; }
+
+ virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
+ { return MB_FAILURE; }
+
+ virtual std::string get_mesh_type_name()
+ { return "CAM_UNKNOWN"; }
+};
+
+//! Child helper class for unknown grid (NOT_CAM)
+class NCHNotCam : public NCHelper
+{
+public:
+ NCHNotCam(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+ { return MB_FAILURE; }
+
+ virtual ErrorCode init_local_gid(EntityHandle file_set)
+ { return MB_FAILURE; }
+
+ virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
+ { return MB_FAILURE; }
+
+ virtual std::string get_mesh_type_name()
+ { return "NOT_CAM"; }
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/ReadNC.cpp b/src/io/ReadNC.cpp
index b194314..42398b9 100644
--- a/src/io/ReadNC.cpp
+++ b/src/io/ReadNC.cpp
@@ -1,4 +1,5 @@
#include "ReadNC.hpp"
+#include "NCHelper.hpp"
#include <algorithm>
#include <assert.h>
@@ -90,6 +91,8 @@ void ReadNC::reset() {
ReadNC::~ReadNC() {
mbImpl->release_interface(readMeshIface);
+ if (helper != NULL)
+ delete helper;
}
ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set, const FileOptions& opts,
@@ -114,7 +117,6 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
mpFileIdTag = file_id_tag; // store the pointer to the tag ; if not null, set when global id tag
// is set too, with the same data , duplicated
-
std::string partition_tag_name;
rval = parse_options(opts, var_names, tstep_nums, tstep_vals);
ERRORR(rval, "Trouble parsing option string.");
@@ -152,7 +154,6 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
// end of BIL
-
// Read the header (num dimensions, dimensions, num variables, global attribs)
rval = read_header();
ERRORR(rval, " ");
@@ -189,8 +190,11 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
delete helper;
helper = NCHelper::get_nc_helper(this, fileId, opts);
+ if (helper == NULL) {
+ ERRORR(MB_FAILURE, "Failed to get NCHelper class instance.");
+ }
- rval = helper->init_vals(opts, tmp_set);
+ rval = helper->init_mesh_vals(opts, tmp_set);
if (MB_SUCCESS != rval)
return rval;
@@ -201,38 +205,21 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
ERRORR(rval, "Mesh characteristics didn't match from last read.\n");
}
else if (!noMesh) {
- if (CAM_SE == helper->get_cam_type())
- rval = create_ucd_verts_quads(opts, tmp_set, quads);
- else
- rval = create_verts_quads(scdi, tmp_set, quads);
+ rval = helper->create_verts_quads(scdi, opts, tmp_set, quads);
ERRORR(rval, "Trouble creating vertices and quads.");
}
- if (noMesh && CAM_SE == helper->get_cam_type())
+
+ if (noMesh)
{
- // we need to populate localGid range with the gids of vertices from the tmp_set
- // localGid is important in reading the variable data into the nodes
- // also, for our purposes, localGid is truly the GLOBAL_ID tag data, not other
- // file_id tags that could get passed around in other scenarios for parallel reading
- // for nodal_partition, this local gid is easier, should be initialized with only
- // the owned nodes
-
- // we need to get all vertices from tmp_set (it is the input set in no_mesh scenario)
- Range local_verts;
- rval = mbImpl->get_entities_by_dimension(tmp_set, 0, local_verts);
- if (MB_FAILURE == rval)
- return rval;
- std::vector<int> gids(local_verts.size());
- // !IMPORTANT : this has to be the GLOBAL_ID tag
- rval=mbImpl->tag_get_data(mGlobalIdTag, local_verts, &gids[0]);
- if (MB_FAILURE == rval)
+ // Initialize local gid for HOMME
+ rval = helper->init_local_gid(tmp_set);
+ if (MB_SUCCESS != rval)
return rval;
- // this will do a smart copy
- std::copy(gids.begin(), gids.end(), range_inserter(localGid));
}
// Read variables onto grid
if (!noVars) {
- rval = read_variables(tmp_set, var_names, tstep_nums);
+ rval = helper->read_variables(tmp_set, var_names, tstep_nums);
if (MB_FAILURE == rval)
return rval;
}
@@ -244,7 +231,7 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
if (mit != varInfo.end())
filteredDimNames.push_back(dimNames[i]);
}
- rval = read_variables(tmp_set, filteredDimNames, tstep_nums);
+ rval = helper->read_variables(tmp_set, filteredDimNames, tstep_nums);
if (MB_FAILURE == rval)
return rval;
}
@@ -297,72 +284,7 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
return MB_SUCCESS;
}
-ErrorCode ReadNC::check_conventions_attribute()
-{
- isCf = false;
-
- std::map<std::string, AttData>::iterator attIt = globalAtts.find("conventions");
- if (attIt == globalAtts.end())
- attIt = globalAtts.find("Conventions");
-
- if (attIt == globalAtts.end()) {
- ERRORR(MB_FAILURE, "File does not have conventions global attribute.\n");
- }
-
- unsigned int sz = attIt->second.attLen;
- std::string att_data;
- att_data.resize(sz + 1);
- att_data[sz] = '\000';
- int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
- ERRORS(success, "Failed to read conventions global attribute char data.");
- if (att_data.find("CF") == std::string::npos) {
- ERRORR(MB_FAILURE, "File not following known conventions.\n");
- }
- else
- isCf = true;
-
- return MB_SUCCESS;
-}
-
-ErrorCode ReadNC::check_source_attribute()
-{
- isCam = false;
-
- std::map<std::string, AttData>::iterator attIt = globalAtts.find("source");
-
- if (attIt == globalAtts.end()) {
- ERRORR(MB_FAILURE, "File does not have source global attribute.\n");
- }
-
- unsigned int sz = attIt->second.attLen;
- std::string att_data;
- att_data.resize(sz + 1);
- att_data[sz] = '\000';
- int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
- ERRORS(success, "Failed to read source global attribute char data.");
- if (att_data.find("CAM") != std::string::npos)
- isCam = true;
-
- return MB_SUCCESS;
-}
-
-ErrorCode ReadNC::check_np_attribute()
-{
- std::map<std::string, AttData>::iterator attIt = globalAtts.find("np");
- if (attIt != globalAtts.end())
- {
- int success = NCFUNC(get_att_int)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &spectralOrder);
- ERRORS(success, "Failed to read np global attribute int data.");
- spectralOrder--; // Spectral order is one less than np
-
- return MB_SUCCESS;
- }
-
- return MB_FAILURE;
-}
-
ErrorCode ReadNC::load_BIL(std::string , const EntityHandle* , const FileOptions& , const Tag* ) {
-
/*
BIL_Init( MPI_COMM_WORLD );
@@ -514,7 +436,6 @@ ErrorCode ReadNC::parse_options(const FileOptions &opts, std::vector<std::string
partMethod = ScdParData::ALLJORKORI;
else
partMethod = dum;
-
#endif
return MB_SUCCESS;
@@ -545,7 +466,7 @@ ErrorCode ReadNC::check_verts_quads(EntityHandle file_set) {
return MB_SUCCESS;
}
-ErrorCode ReadNC::create_verts_quads(ScdInterface *scdi, EntityHandle tmp_set, Range &quads) {
+ErrorCode ReadNC::create_scd_verts_quads(ScdInterface *scdi, EntityHandle file_set, Range &quads) {
Range tmp_range;
ScdBox *scd_box;
@@ -557,7 +478,7 @@ ErrorCode ReadNC::create_verts_quads(ScdInterface *scdi, EntityHandle tmp_set, R
tmp_range.insert(scd_box->start_vertex(), scd_box->start_vertex() + scd_box->num_vertices() - 1);
tmp_range.insert(scd_box->start_element(), scd_box->start_element() + scd_box->num_elements() - 1);
tmp_range.insert(scd_box->box_set());
- rval = mbImpl->add_entities(tmp_set, tmp_range);
+ rval = mbImpl->add_entities(file_set, tmp_range);
ERRORR(rval, "Couldn't add new vertices to file set.");
dbgOut.tprintf(1, "scdbox %d quads, %d vertices\n", scd_box->num_elements(), scd_box->num_vertices());
@@ -636,7 +557,7 @@ ErrorCode ReadNC::create_verts_quads(ScdInterface *scdi, EntityHandle tmp_set, R
return MB_SUCCESS;
}
-ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle tmp_set, Range &quads) {
+ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle file_set, Range &quads) {
// need to get/read connectivity data before creating elements
std::string conn_fname;
@@ -733,7 +654,7 @@ ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle t
// start_idx is the starting index in the HommeMapping connectivity list for this proc, before converting to coarse quad representation
start_idx = 4 * rank * num_coarse_quads * spectral_unit;
// iextra = # coarse quads extra after equal split over procs
- int iextra = num_quads % (procs*spectral_unit);
+ int iextra = num_quads % (procs*spectral_unit);
if (rank < iextra) num_coarse_quads++;
start_idx += 4 * spectral_unit * std::min(rank, iextra);
// num_fine_quads is the number of quads in the connectivity list in HommeMapping file assigned to this proc
@@ -749,7 +670,7 @@ ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle t
SpectralMeshTool smt(mbImpl, spectralOrder);
if (!spectralMesh) {
rval = readMeshIface->get_element_connect(num_coarse_quads, 4,
- MBQUAD, 0, start_quad, conn_arr,
+ MBQUAD, 0, start_quad, conn_arr,
// might have to create gather mesh later
(create_gathers ? num_coarse_quads + num_quads : num_coarse_quads));
ERRORR(rval, "Failed to create quads.");
@@ -763,18 +684,18 @@ ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle t
int count, v_per_e;
rval = mbImpl->connect_iterate(tmp_range.begin(), tmp_range.end(), conn_arr, v_per_e, count);
ERRORR(rval, "Failed to get connectivity of spectral elements.");
- rval = mbImpl->tag_iterate(smt.spectral_vertices_tag(true), tmp_range.begin(), tmp_range.end(),
+ rval = mbImpl->tag_iterate(smt.spectral_vertices_tag(true), tmp_range.begin(), tmp_range.end(),
count, (void*&)sv_ptr);
ERRORR(rval, "Failed to get fine connectivity of spectral elements.");
}
-
+
// on this proc, I get columns ldims[1]..ldims[4], inclusive; need to find which vertices those correpond to
unsigned int num_local_verts = localGid.size();
unsigned int num_total_verts = gDims[3] - gDims[0] + 1;
// create vertices
std::vector<double*> arrays;
- rval = readMeshIface->get_node_coords(3, num_local_verts, 0, start_vertex, arrays,
+ rval = readMeshIface->get_node_coords(3, num_local_verts, 0, start_vertex, arrays,
// might have to create gather mesh later
(create_gathers ? num_local_verts+num_total_verts : num_local_verts));
ERRORR(rval, "Couldn't create vertices in ucd mesh.");
@@ -840,14 +761,14 @@ ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle t
// add new vertices and elements to the set
quads.merge(tmp_range);
tmp_range.insert(start_vertex, start_vertex + num_local_verts - 1);
- rval = mbImpl->add_entities(tmp_set, tmp_range);
+ rval = mbImpl->add_entities(file_set, tmp_range);
ERRORR(rval, "Couldn't add new vertices and quads/hexes to file set.");
// mark the set with the spectral order
Tag sporder;
rval = mbImpl->tag_get_handle("SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, sporder, MB_TAG_CREAT | MB_TAG_SPARSE);
ERRORR(rval, "Couldn't create spectral order tag.");
- rval = mbImpl->tag_set_data(sporder, &tmp_set, 1, &spectralOrder);
+ rval = mbImpl->tag_set_data(sporder, &file_set, 1, &spectralOrder);
ERRORR(rval, "Couldn't set value for spectral order tag.");
/*
bool gatherOpt = false;
@@ -855,7 +776,7 @@ ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle t
gatherOpt = true;
*/
/*
- - if (root)
+ - if (root)
. create vertices for all vertex positions in 2D grid
. create quads for all vertex positions
. create new entity set & insert these vertices/quads into it
@@ -866,16 +787,16 @@ ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle t
#ifdef USE_MPI
if (isParallel && myPcomm->proc_config().proc_rank() == 0) {
#endif
- EntityHandle gather_set;
+ EntityHandle gather_set;
rval = mbImpl->create_meshset(MESHSET_SET, gather_set);
ERRORR(rval, "Trouble creating gather set.");
// create vertices
arrays.clear();
// don't need to specify allocation number here, because we know enough verts were created before
- rval = readMeshIface->get_node_coords(3, num_total_verts, 0, start_vertex, arrays);
+ rval = readMeshIface->get_node_coords(3, num_total_verts, 0, start_vertex, arrays);
ERRORR(rval, "Couldn't create vertices in ucd mesh for gather set.");
-
+
xptr = arrays[0], yptr = arrays[1], zptr = arrays[2];
for (i = 0; i < (int)num_total_verts; ++i) {
double cosphi = cos(pideg * jlVals[i]);
@@ -887,9 +808,9 @@ ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle t
yptr[i] = rad * ymult;
zptr[i] = rad * zmult;
}
-
+
// get ptr to gid memory for vertices
- Range gather_verts(start_vertex, start_vertex + num_total_verts - 1);
+ Range gather_verts(start_vertex, start_vertex + num_total_verts - 1);
rval = mbImpl->tag_iterate(mGlobalIdTag, gather_verts.begin(), gather_verts.end(), count, data);
ERRORR(rval, "Failed to get tag iterator.");
assert(count == (int) num_total_verts);
@@ -923,13 +844,13 @@ ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle t
ERRORR(rval, "Couldn't add quads to gather set.");
Tag gathersettag;
- rval = mbImpl->tag_get_handle("GATHER_SET", 1, MB_TYPE_INTEGER, gathersettag,
+ rval = mbImpl->tag_get_handle("GATHER_SET", 1, MB_TYPE_INTEGER, gathersettag,
MB_TAG_CREAT | MB_TAG_SPARSE);
ERRORR(rval, "Couldn't create gather set tag.");
int gatherval = 1;
rval = mbImpl->tag_set_data(gathersettag, &gather_set, 1, &gatherval);
ERRORR(rval, "Couldn't set value for gather set tag.");
-
+
#ifdef USE_MPI
}
#endif
@@ -938,10 +859,10 @@ ErrorCode ReadNC::create_ucd_verts_quads(const FileOptions &opts, EntityHandle t
}
ErrorCode ReadNC::read_variable_setup(std::vector<std::string> &var_names, std::vector<int> &tstep_nums,
- std::vector<VarData> &vdatas, std::vector<VarData> &vsetdatas) {
+ std::vector<VarData> &vdatas, std::vector<VarData> &vsetdatas, bool is_scd_mesh) {
std::map<std::string, VarData>::iterator mit;
- if (helper->get_cam_type() != CAM_SE) { // scd mesh
+ if (is_scd_mesh) { // scd mesh
// if empty read them all
if (var_names.empty()) {
for (mit = varInfo.begin(); mit != varInfo.end(); mit++) {
@@ -1042,7 +963,7 @@ ErrorCode ReadNC::read_variable_setup(std::vector<std::string> &var_names, std::
return MB_SUCCESS;
}
-ErrorCode ReadNC::read_variable_allocate(EntityHandle file_set, std::vector<VarData> &vdatas, std::vector<int> &tstep_nums) {
+ErrorCode ReadNC::read_variable_allocate(EntityHandle file_set, std::vector<VarData> &vdatas, std::vector<int> &tstep_nums, bool is_scd_mesh) {
ErrorCode rval = MB_SUCCESS;
std::vector<EntityHandle>* ehandles = NULL;
@@ -1140,7 +1061,7 @@ ErrorCode ReadNC::read_variable_allocate(EntityHandle file_set, std::vector<VarD
switch (vdatas[i].entLoc) {
case 0:
// vertices
- if (helper->get_cam_type() != CAM_SE) {
+ if (is_scd_mesh) {
// only structured mesh has j parameter that multiplies i to get total # vertices
vdatas[i].readDims[t].push_back(lDims[1]);
vdatas[i].readCounts[t].push_back(lDims[4] - lDims[1] + 1);
@@ -1242,33 +1163,35 @@ ErrorCode ReadNC::read_variable_allocate(EntityHandle file_set, std::vector<VarD
return rval;
}
-ErrorCode ReadNC::read_variables(EntityHandle file_set, std::vector<std::string> &var_names, std::vector<int> &tstep_nums) {
+ErrorCode ReadNC::read_variables(EntityHandle file_set, std::vector<std::string> &var_names, std::vector<int> &tstep_nums, bool is_scd_mesh) {
std::vector<VarData> vdatas;
std::vector<VarData> vsetdatas;
- ErrorCode rval = read_variable_setup(var_names, tstep_nums, vdatas, vsetdatas);
+ ErrorCode rval = read_variable_setup(var_names, tstep_nums, vdatas, vsetdatas, is_scd_mesh);
ERRORR(rval, "Trouble setting up read variable.");
- if (helper->get_cam_type() != CAM_SE) {
+ if (is_scd_mesh) {
// create COORDS tag for quads
rval = create_quad_coordinate_tag(file_set);
ERRORR(rval, "Trouble creating coordinate tags to entities quads");
}
if (!vsetdatas.empty()) {
- rval = read_variable_to_set(file_set, vsetdatas, tstep_nums);
+ rval = read_variable_to_set(file_set, vsetdatas, tstep_nums, is_scd_mesh);
ERRORR(rval, "Trouble read variables to set.");
}
if (!vdatas.empty()) {
#ifdef PNETCDF_FILE
- if (helper->get_cam_type() == CAM_SE) // in serial, we will use the old read, everything is contiguous
+ if (!is_scd_mesh) // in serial, we will use the old read, everything is contiguous
// in parallel, we will use async read in pnetcdf
// the other mechanism is not working, forget about it
rval = read_variable_to_nonset_async(file_set, vdatas, tstep_nums);
else
+ rval = read_variable_to_nonset(file_set, vdatas, tstep_nums, is_scd_mesh);
+#else
+ rval = read_variable_to_nonset(file_set, vdatas, tstep_nums, is_scd_mesh);
#endif
- rval = read_variable_to_nonset(file_set, vdatas, tstep_nums);
ERRORR(rval, "Trouble read variables to entities verts/edges/quads.");
}
@@ -1276,7 +1199,6 @@ ErrorCode ReadNC::read_variables(EntityHandle file_set, std::vector<std::string>
return MB_SUCCESS;
}
-
ErrorCode ReadNC::read_variable_to_set_allocate(std::vector<VarData> &vdatas, std::vector<int> &tstep_nums) {
ErrorCode rval = MB_SUCCESS;
@@ -1357,7 +1279,7 @@ ErrorCode ReadNC::read_variable_to_set_allocate(std::vector<VarData> &vdatas, st
return rval;
}
-ErrorCode ReadNC::read_variable_to_set(EntityHandle file_set, std::vector<VarData> &vdatas, std::vector<int> &tstep_nums) {
+ErrorCode ReadNC::read_variable_to_set(EntityHandle file_set, std::vector<VarData> &vdatas, std::vector<int> &tstep_nums, bool is_scd_mesh) {
ErrorCode rval = read_variable_to_set_allocate(vdatas, tstep_nums);
ERRORR(rval, "Trouble allocating read variables to set.");
@@ -1421,7 +1343,7 @@ ErrorCode ReadNC::read_variable_to_set(EntityHandle file_set, std::vector<VarDat
for (unsigned int i = 0; i < vdatas.size(); i++) {
for (unsigned int t = 0; t < tstep_nums.size(); t++) {
dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
- ErrorCode tmp_rval = convert_variable(vdatas[i], t);
+ ErrorCode tmp_rval = convert_variable(vdatas[i], t, is_scd_mesh);
if (MB_SUCCESS != tmp_rval)
rval = tmp_rval;
if (vdatas[i].varDims.size() <= 1)
@@ -1450,8 +1372,8 @@ ErrorCode ReadNC::read_variable_to_set(EntityHandle file_set, std::vector<VarDat
return rval;
}
-ErrorCode ReadNC::read_variable_to_nonset(EntityHandle file_set, std::vector<VarData> &vdatas, std::vector<int> &tstep_nums) {
- ErrorCode rval = read_variable_allocate(file_set, vdatas, tstep_nums);
+ErrorCode ReadNC::read_variable_to_nonset(EntityHandle file_set, std::vector<VarData> &vdatas, std::vector<int> &tstep_nums, bool is_scd_mesh) {
+ ErrorCode rval = read_variable_allocate(file_set, vdatas, tstep_nums, is_scd_mesh);
ERRORR(rval, "Trouble allocating read variables.");
// finally, read into that space
@@ -1462,7 +1384,7 @@ ErrorCode ReadNC::read_variable_to_nonset(EntityHandle file_set, std::vector<Var
void *data = vdatas[i].varDatas[t];
std::size_t sz = 1;
size_t ni = vdatas[i].readCounts[t][2], nj = vdatas[i].readCounts[t][3], nk = vdatas[i].readCounts[t][1];
- if (helper->get_cam_type() != CAM_SE) {
+ if (is_scd_mesh) {
for (std::size_t idx = 0; idx != vdatas[i].readCounts[t].size(); ++idx)
sz *= vdatas[i].readCounts[t][idx];
}
@@ -1504,7 +1426,7 @@ ErrorCode ReadNC::read_variable_to_nonset(EntityHandle file_set, std::vector<Var
case NC_FLOAT: {
std::vector<float> tmpfloatdata(sz);
- if (helper->get_cam_type() != CAM_SE)
+ if (is_scd_mesh)
{
success = NCFUNCAG(_vara_float)(fileId, vdatas[i].varId, &vdatas[i].readDims[t][0], &vdatas[i].readCounts[t][0],
&tmpfloatdata[0] NCREQ);
@@ -1600,7 +1522,7 @@ ErrorCode ReadNC::read_variable_to_nonset(EntityHandle file_set, std::vector<Var
for (unsigned int i = 0; i < vdatas.size(); i++) {
for (unsigned int t = 0; t < tstep_nums.size(); t++) {
dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
- ErrorCode tmp_rval = convert_variable(vdatas[i], t);
+ ErrorCode tmp_rval = convert_variable(vdatas[i], t, is_scd_mesh);
if (MB_SUCCESS != tmp_rval)
rval = tmp_rval;
}
@@ -1621,7 +1543,7 @@ ErrorCode ReadNC::read_variable_to_nonset(EntityHandle file_set, std::vector<Var
ErrorCode ReadNC::read_variable_to_nonset_async(EntityHandle file_set, std::vector<VarData> &vdatas,
std::vector<int> &tstep_nums)
{
- ErrorCode rval = read_variable_allocate(file_set, vdatas, tstep_nums);
+ ErrorCode rval = read_variable_allocate(file_set, vdatas, tstep_nums, false);
ERRORR(rval, "Trouble allocating read variables.");
// finally, read into that space
@@ -1692,7 +1614,6 @@ ErrorCode ReadNC::read_variable_to_nonset_async(EntityHandle file_set, std::vect
success = ncmpi_wait_all(fileId, requests.size(), &requests[0], &statuss[0]);
ERRORS(success, "Failed on wait_all.");
-
if (vdatas[i].numLev != 1)
// switch from k varying slowest to k varying fastest
success = kji_to_jik_stride(ni, nj, nk, data, &tmpdoubledata[0]);
@@ -1773,7 +1694,7 @@ ErrorCode ReadNC::read_variable_to_nonset_async(EntityHandle file_set, std::vect
for (unsigned int i = 0; i < vdatas.size(); i++) {
for (unsigned int t = 0; t < tstep_nums.size(); t++) {
dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
- ErrorCode tmp_rval = convert_variable(vdatas[i], t);
+ ErrorCode tmp_rval = convert_variable(vdatas[i], t, false);
if (MB_SUCCESS != tmp_rval)
rval = tmp_rval;
}
@@ -1790,12 +1711,12 @@ ErrorCode ReadNC::read_variable_to_nonset_async(EntityHandle file_set, std::vect
}
#endif
-ErrorCode ReadNC::convert_variable(VarData &var_data, int tstep_num) {
+ErrorCode ReadNC::convert_variable(VarData &var_data, int tstep_num, bool is_scd_mesh) {
// get ptr to tag space
void *data = var_data.varDatas[tstep_num];
std::size_t sz = 1;
- if (helper->get_cam_type() != CAM_SE) {
+ if (is_scd_mesh) {
for (std::size_t idx = 0; idx != var_data.readCounts[tstep_num].size(); ++idx)
sz *= var_data.readCounts[tstep_num][idx];
}
@@ -2364,7 +2285,7 @@ ErrorCode ReadNC::init_EulSpcscd_vals(const FileOptions &opts, EntityHandle file
gDims[2] = -1;
gDims[5] = -1;
- // look for time dimensions
+ // look for time dimensions
if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
idx = vit - dimNames.begin();
else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
@@ -3250,17 +3171,7 @@ ErrorCode ReadNC::create_tags(ScdInterface *scdi, EntityHandle file_set, const s
// <__MESH_TYPE>
Tag meshTypeTag = 0;
tag_name = "__MESH_TYPE";
- std::string meshTypeName;
- switch (helper->get_cam_type())
- {
- case CAM_EUL: meshTypeName="CAM_EUL"; break;
- case CAM_FV: meshTypeName="CAM_FV"; break;
- case CAM_SE: meshTypeName="CAM_SE"; break;
- case CAM_UNKNOWN: meshTypeName="CAM_UNKNOWN"; break;
- case NOT_CAM: meshTypeName="NOT_CAM"; break;
- default: meshTypeName="NOT_CAM"; break;
-
- }
+ std::string meshTypeName = helper->get_mesh_type_name();
rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, meshTypeTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
ERRORR(rval, "Trouble creating __MESH_TYPE tag.");
@@ -3395,107 +3306,93 @@ ErrorCode ReadNC::create_quad_coordinate_tag(EntityHandle file_set) {
return MB_SUCCESS;
}
-NCHelper* NCHelper::get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts)
+ErrorCode ReadNC::check_conventions_attribute()
{
- // Unknown grid, will fill this in later for POP, CICE and CLM
- if (!readNC->isCam)
- return new NCHNotCam(readNC, fileId);
-
- // If a CAM file, which type?
- if (NCHEuler::can_read_file(readNC))
- return new NCHEuler(readNC, fileId);
- else if (NCHFV::can_read_file(readNC))
- return new NCHFV(readNC, fileId);
- else if (NCHHomme::can_read_file(readNC, opts))
- return new NCHHomme(readNC, fileId);
-
- // Unknown CAM grid
- return new NCHUnknownCam(readNC, fileId);
-}
+ isCf = false;
-bool NCHEuler::can_read_file(ReadNC* readNC)
-{
- std::vector<std::string>& dimNames = readNC->dimNames;
+ std::map<std::string, AttData>::iterator attIt = globalAtts.find("conventions");
+ if (attIt == globalAtts.end())
+ attIt = globalAtts.find("Conventions");
- // If dimension names "lon" AND "lat' exist then it's the Eulerian Spectral grid or the FV grid
- if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
- dimNames.end(), std::string("lat")) != dimNames.end()))
- {
- // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
- if ((std::find(dimNames.begin(), dimNames.end(), std::string("slon")) != dimNames.end()) && (std::find(dimNames.begin(),
- dimNames.end(), std::string("slat")) != dimNames.end()))
- return false;
- else
- return true;
+ if (attIt == globalAtts.end()) {
+ ERRORR(MB_FAILURE, "File does not have conventions global attribute.\n");
}
- return false;
-}
-
-ErrorCode NCHEuler::init_vals(const FileOptions& opts, EntityHandle file_set)
-{
- ErrorCode rval = _readNC->init_EulSpcscd_vals(opts, file_set);
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble initializing Euler grid.");
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ ERRORS(success, "Failed to read conventions global attribute char data.");
+ if (att_data.find("CF") == std::string::npos) {
+ ERRORR(MB_FAILURE, "File not following known conventions.\n");
+ }
+ else
+ isCf = true;
- return rval;
+ return MB_SUCCESS;
}
-bool NCHFV::can_read_file(ReadNC* readNC)
+ErrorCode ReadNC::check_source_attribute()
{
- std::vector<std::string>& dimNames = readNC->dimNames;
+ isCam = false;
- // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
- if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
- dimNames.end(), std::string("lat")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slon"))
- != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end()))
- return true;
+ std::map<std::string, AttData>::iterator attIt = globalAtts.find("source");
- return false;
-}
+ if (attIt == globalAtts.end()) {
+ ERRORR(MB_FAILURE, "File does not have source global attribute.\n");
+ }
-ErrorCode NCHFV::init_vals(const FileOptions& opts, EntityHandle file_set)
-{
- ErrorCode rval = _readNC->init_FVCDscd_vals(opts, file_set);
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble initializing FV grid.");
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ ERRORS(success, "Failed to read source global attribute char data.");
+ if (att_data.find("CAM") != std::string::npos)
+ isCam = true;
- return rval;
+ return MB_SUCCESS;
}
-bool NCHHomme::can_read_file(ReadNC* readNC, const FileOptions& opts)
+ErrorCode ReadNC::init_local_gid(EntityHandle file_set)
{
- // If global attribute "np" exists then it's the HOMME grid
- ErrorCode rval = readNC->check_np_attribute();
- if (MB_SUCCESS == rval) {
- if (MB_SUCCESS == opts.match_option("PARTITION_METHOD", "NODAL_PARTITION"))
- readNC->partMethod = -1;
-
- return true;
- }
+ // we need to populate localGid range with the gids of vertices from the file_set
+ // localGid is important in reading the variable data into the nodes
+ // also, for our purposes, localGid is truly the GLOBAL_ID tag data, not other
+ // file_id tags that could get passed around in other scenarios for parallel reading
+ // for nodal_partition, this local gid is easier, should be initialized with only
+ // the owned nodes
+
+ // we need to get all vertices from file_set (it is the input set in no_mesh scenario)
+ Range local_verts;
+ ErrorCode rval = mbImpl->get_entities_by_dimension(file_set, 0, local_verts);
+ if (MB_FAILURE == rval)
+ return rval;
- return false;
-}
+ std::vector<int> gids(local_verts.size());
+ // !IMPORTANT : this has to be the GLOBAL_ID tag
+ rval = mbImpl->tag_get_data(mGlobalIdTag, local_verts, &gids[0]);
+ if (MB_FAILURE == rval)
+ return rval;
-ErrorCode NCHHomme::init_vals(const FileOptions& opts, EntityHandle file_set)
-{
- ErrorCode rval = _readNC->init_HOMMEucd_vals();
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble initializing Homme grid.");
+ // this will do a smart copy
+ std::copy(gids.begin(), gids.end(), range_inserter(localGid));
- return rval;
+ return MB_SUCCESS;
}
-ErrorCode NCHUnknownCam::init_vals(const FileOptions& opts, EntityHandle file_set)
+ErrorCode ReadNC::check_np_attribute()
{
- _readNC->readMeshIface->report_error("%s", "Unknown CAM grid.");
-
- return MB_FAILURE;
-}
+ std::map<std::string, AttData>::iterator attIt = globalAtts.find("np");
+ if (attIt != globalAtts.end())
+ {
+ int success = NCFUNC(get_att_int)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &spectralOrder);
+ ERRORS(success, "Failed to read np global attribute int data.");
+ spectralOrder--; // Spectral order is one less than np
-ErrorCode NCHNotCam::init_vals(const FileOptions& opts, EntityHandle file_set)
-{
- _readNC->readMeshIface->report_error("%s", "Unknown grid.");
+ return MB_SUCCESS;
+ }
return MB_FAILURE;
}
diff --git a/src/io/ReadNC.hpp b/src/io/ReadNC.hpp
index cf6fb7d..b8d45fd 100644
--- a/src/io/ReadNC.hpp
+++ b/src/io/ReadNC.hpp
@@ -64,7 +64,6 @@ class ReadUtilIface;
class ScdInterface;
class NCHelper;
-
//! Output Exodus File for VERDE
class ReadNC : public ReaderIface
{
@@ -76,9 +75,9 @@ class ReadNC : public ReaderIface
friend class NCHNotCam;
public:
-
+
static ReaderIface* factory( Interface* );
-
+
//! load an NC file
ErrorCode load_file( const char* file_name,
const EntityHandle* file_set,
@@ -86,8 +85,8 @@ public:
const SubsetList* subset_list = 0,
const Tag* file_id_tag = 0 );
- //! Constructor
- ReadNC(Interface* impl = NULL);
+ //! Constructor
+ ReadNC(Interface* impl = NULL);
//! Destructor
virtual ~ReadNC();
@@ -97,12 +96,10 @@ public:
const FileOptions& opts,
std::vector<int>& tag_values_out,
const SubsetList* subset_list = 0 );
-
+
// ENTLOCNSEDGE for north/south edge
// ENTLOCWEEDGE for west/east edge
- enum EntityLocation {ENTLOCNODE=0, ENTLOCNSEDGE, ENTLOCEWEDGE, ENTLOCQUAD, ENTLOCSET};
-
- enum CamType {CAM_EUL=0, CAM_FV, CAM_SE, CAM_UNKNOWN, NOT_CAM};
+ enum EntityLocation {ENTLOCNODE = 0, ENTLOCNSEDGE, ENTLOCEWEDGE, ENTLOCQUAD, ENTLOCSET};
private:
@@ -141,7 +138,7 @@ private:
ReadUtilIface* readMeshIface;
bool dimension_exists(const char *attrib_name);
-
+
void reset();
//! read the header information
@@ -150,28 +147,25 @@ private:
//! get all global attributes in the file
ErrorCode get_attributes(int var_id, int num_atts, std::map<std::string,AttData> &atts,
const char *prefix="");
-
+
//! get all dimensions in the file
ErrorCode get_dimensions(int file_id, std::vector<std::string> &dim_names, std::vector<int> &dim_vals);
//! get the variable names and other info defined for this file
ErrorCode get_variables();
-
+
//! parse min/max i/j/k in options, if any
ErrorCode init_EulSpcscd_vals(const FileOptions &opts, EntityHandle file_set);
ErrorCode init_FVCDscd_vals(const FileOptions &opts, EntityHandle file_set);
ErrorCode read_coordinate(const char *var_name, int lmin, int lmax,
std::vector<double> &cvals);
-
- ErrorCode read_coordinate_nc(const char *var_name, int lmin, int lmax,
- std::vector<double> &cvals);
-
+
//! number of dimensions in this nc file
unsigned int number_dimensions();
//! create vertices for the file
- ErrorCode create_verts_quads(ScdInterface *scdi, EntityHandle file_set, Range &quads);
+ ErrorCode create_scd_verts_quads(ScdInterface *scdi, EntityHandle file_set, Range &quads);
//! check number of vertices and elements against what's already in file_set
ErrorCode check_verts_quads(EntityHandle file_set);
@@ -183,54 +177,40 @@ private:
ErrorCode read_variable_to_set_allocate(std::vector<VarData> &vdatas,
std::vector<int> &tstep_nums);
-
+
ErrorCode read_variable_to_set(EntityHandle file_set, std::vector<VarData> &vdatas,
- std::vector<int> &tstep_nums);
-
+ std::vector<int> &tstep_nums, bool is_scd_mesh);
+
ErrorCode read_variable_to_nonset(EntityHandle file_set, std::vector<VarData> &vdatas,
- std::vector<int> &tstep_nums);
+ std::vector<int> &tstep_nums, bool is_scd_mesh);
#ifdef PNETCDF_FILE
ErrorCode read_variable_to_nonset_async(EntityHandle file_set, std::vector<VarData> &vdatas,
std::vector<int> &tstep_nums);
#endif
- ErrorCode read_variable_to_ucdmesh(EntityHandle file_set,
- std::vector<std::string> &varnames,
- std::vector<int> &tstep_nums,
- std::vector<VarData> &vdatas);
-
- ErrorCode read_variable_ucd_allocate(std::vector<VarData> &vdatas,
- std::vector<int> &tstep_nums,
- Range &verts);
-
- ErrorCode get_tag_ucd(VarData &var_data, int tstep_num, Tag &tagh);
-
- ErrorCode read_variable_ucd_setup(std::vector<std::string> &var_names,
- std::vector<int> &tstep_nums,
- std::vector<VarData> &vdatas);
-
ErrorCode read_variables(EntityHandle file_set, std::vector<std::string> &var_names,
- std::vector<int> &tstep_nums);
-
+ std::vector<int> &tstep_nums, bool is_scd_mesh);
+
ErrorCode read_variable_allocate(EntityHandle file_set, std::vector<VarData> &vdatas,
- std::vector<int> &tstep_nums);
-
+ std::vector<int> &tstep_nums, bool is_scd_mesh);
+
ErrorCode read_variable_setup(std::vector<std::string> &var_names,
std::vector<int> &tstep_nums,
std::vector<VarData> &vdatas,
- std::vector<VarData> &vsetdatas);
-
- ErrorCode convert_variable(VarData &var_data, int tstep_num);
-
+ std::vector<VarData> &vsetdatas,
+ bool is_scd_mesh);
+
+ ErrorCode convert_variable(VarData &var_data, int tstep_num, bool is_scd_mesh);
+
ErrorCode get_tag_to_set(VarData &var_data, int tstep_num, Tag &tagh);
-
+
ErrorCode get_tag(VarData &var_data, int tstep_num, Tag &tagh, int num_lev);
-
+
//! create nc conventional tags
ErrorCode create_tags(ScdInterface *scdi, EntityHandle file_set,
const std::vector<int> &tstep_nums);
-
+
//! create a character string attString of attMap. with '\0'
//! terminating each attribute name, ';' separating the data type
//! and value, and ';' separating one name/data type/value from
@@ -249,8 +229,8 @@ private:
ErrorCode init_HOMMEucd_vals();
- ErrorCode create_ucd_verts_quads(const FileOptions &opts, EntityHandle tmp_set, Range &quads);
-
+ ErrorCode create_ucd_verts_quads(const FileOptions &opts, EntityHandle file_set, Range &quads);
+
ErrorCode load_BIL(std::string dir_name,
const EntityHandle* file_set,
const FileOptions& opts,
@@ -260,12 +240,6 @@ private:
bool BIL_mode_enabled(const char * file_name);
- ErrorCode check_conventions_attribute();
-
- ErrorCode check_source_attribute();
-
- ErrorCode check_np_attribute();
-
template <typename T> ErrorCode kji_to_jik(size_t ni, size_t nj, size_t nk, void *dest, T *source)
{
size_t nik = ni * nk, nij = ni * nj;
@@ -276,13 +250,12 @@ private:
tmp_data[j*nik+i*nk+k] = source[k*nij+j*ni+i];
return MB_SUCCESS;
}
-
+
// this version takes as input the moab range, from which we actually need just the
// size of each sequence, for a proper transpose of the data
// we read one time step, one variable at a time, usually, so we will
template <typename T> ErrorCode kji_to_jik_stride(size_t , size_t nj, size_t nk, void *dest, T *source)
{
-
std::size_t idxInSource=0;// position of the start of the stride
// for each subrange, we will transpose a matrix of size subrange*nj*nk (subrange takes
// the role of ni)
@@ -302,11 +275,19 @@ private:
}
return MB_SUCCESS;
}
+
+ ErrorCode check_conventions_attribute();
+
+ ErrorCode check_source_attribute();
+
+ ErrorCode check_np_attribute();
+
+ ErrorCode init_local_gid(EntityHandle file_set);
//------------member variables ------------//
//! interface instance
Interface* mbImpl;
-
+
int CPU_WORD_SIZE;
int IO_WORD_SIZE;
@@ -315,7 +296,7 @@ private:
//! file numbers assigned by netcdf
int fileId, connectId;
-
+
//! dimensions
std::vector<std::string> dimNames;
// these should be taken out when we fix the dummy var info things
@@ -326,13 +307,13 @@ private:
//! global attribs
std::map<std::string,AttData> globalAtts;
-
+
//! variable info
std::map<std::string,VarData> varInfo;
-
+
//! dimensions of grid in file
int gDims[6], tMin, tMax;
-
+
//! dimensions of my part of grid
int lDims[6];
@@ -353,7 +334,7 @@ private:
//! center dimension numbers for i, j
int iCDim, jCDim;
-
+
//! number of the dimension of unlimited dimension, if any
int numUnLim;
@@ -390,14 +371,14 @@ private:
//! partitioning method
int partMethod;
- //! true if a CAM file
+ //! true if a CAM file
bool isCam;
//! true if file is marked as following CF metadata conventions
bool isCf;
-
+
int spectralOrder; // read from variable 'np'
- Range localGid;// used only by camType=CAM_SE
+ Range localGid; // used only by camType = CAM_SE
//! whether mesh is locally periodic in i or j
int locallyPeriodic[2];
@@ -410,23 +391,24 @@ private:
//! directory where data is stored for BIL reader
std::string BIL_dir;
-
+
#ifdef USE_MPI
ParallelComm *myPcomm;
#endif
// read option
bool noMesh;
-
+
// read option
bool noVars;
-
+
// read option
bool spectralMesh;
-
+
// read option
std::string partitionTagName;
+ //! Helper class instance
NCHelper* helper;
};
@@ -436,87 +418,6 @@ inline unsigned int ReadNC::number_dimensions()
return dimVals.size();
}
-// Helper class to isolate reading of several different nc file formats
-class NCHelper
-{
-public:
- NCHelper(ReadNC* readNC, int fileId) : _readNC(readNC), _fileId(fileId) {}
-
- static NCHelper* get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts);
-
- virtual ReadNC::CamType get_cam_type() = 0;
- virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set) = 0;
-
-protected:
- ReadNC* _readNC;
- int _fileId;
-};
-
-// Child helper class for Eulerian Spectral grid (CAM_EUL)
-class NCHEuler : public NCHelper
-{
-public:
- NCHEuler(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
- static bool can_read_file(ReadNC* readNC);
-
-private:
- virtual ReadNC::CamType get_cam_type() { return ReadNC::CAM_EUL; }
- virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
-};
-
-// Child helper class for Finite Volume grid (CAM_FV)
-class NCHFV : public NCHelper
-{
-public:
- NCHFV(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
- static bool can_read_file(ReadNC* readNC);
-
-private:
- virtual ReadNC::CamType get_cam_type() { return ReadNC::CAM_FV; }
- virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
-};
-
-// Child helper class for HOMME grid (CAM_SE)
-class NCHHomme : public NCHelper
-{
-public:
- NCHHomme(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
- static bool can_read_file(ReadNC* readNC, const FileOptions& opts);
-
-private:
- virtual ReadNC::CamType get_cam_type() { return ReadNC::CAM_SE; }
- virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
-};
-
-// Child helper class for unknown CAM grid (CAM_UNKNOWN)
-class NCHUnknownCam : public NCHelper
-{
-public:
- NCHUnknownCam(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
-private:
- virtual ReadNC::CamType get_cam_type() { return ReadNC::CAM_UNKNOWN; }
- virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
-};
-
-// Child helper class for unkown grid (NOT_CAM)
-class NCHNotCam : public NCHelper
-{
-public:
- NCHNotCam(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
-private:
- virtual ReadNC::CamType get_cam_type() { return ReadNC::NOT_CAM; }
- virtual ErrorCode init_vals(const FileOptions& opts, EntityHandle file_set);
-};
-
} // namespace moab
#endif
-
-
-
-
https://bitbucket.org/fathomteam/moab/commits/b749a68a1ff8/
Changeset: b749a68a1ff8
Branch: None
User: danwu
Date: 2013-06-03 15:54:43
Summary: Merged fathomteam/moab into master
Affected #: 4 files
diff --git a/MeshFiles/unittest/io/fv26x46x72.t.3.nc b/MeshFiles/unittest/io/fv26x46x72.t.3.nc
new file mode 100644
index 0000000..15ba5b4
Binary files /dev/null and b/MeshFiles/unittest/io/fv26x46x72.t.3.nc differ
diff --git a/configure.ac b/configure.ac
index c9d6387..130e4a9 100644
--- a/configure.ac
+++ b/configure.ac
@@ -41,6 +41,7 @@ else
AM_SILENT_RULES(no)
fi
])
+AC_C_BIGENDIAN
# Check if platform is BlueGene
AC_MSG_CHECKING([if platform is IBM BlueGene])
diff --git a/src/io/Tqdcfr.cpp b/src/io/Tqdcfr.cpp
index ef01b16..833b15b 100644
--- a/src/io/Tqdcfr.cpp
+++ b/src/io/Tqdcfr.cpp
@@ -39,6 +39,9 @@
#include <sstream>
#include <assert.h>
#include <string.h>
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
namespace moab {
@@ -145,14 +148,54 @@ void Tqdcfr::FREADC( unsigned num_ents ) {
FREADCA( num_ents, &char_buf[0] );
}
+// used for swapping
+static void swap8_voff(long *data)
+{
+ unsigned char tmp, *cdat = (unsigned char *) data;
+ tmp = cdat[0]; cdat[0] = cdat[7], cdat[7] = tmp;
+ tmp = cdat[1]; cdat[1] = cdat[6], cdat[6] = tmp;
+ tmp = cdat[2]; cdat[2] = cdat[5], cdat[5] = tmp;
+ tmp = cdat[3]; cdat[3] = cdat[4], cdat[4] = tmp;
+}
+static void swap4_uint(unsigned int *data)
+{
+ unsigned char tmp, *cdat = (unsigned char *) data;
+ tmp = cdat[0]; cdat[0] = cdat[3], cdat[3] = tmp;
+ tmp = cdat[1]; cdat[1] = cdat[2], cdat[2] = tmp;
+}
+/*static void swap2_ushort(unsigned short *data)
+{
+ unsigned char tmp, *cdat = (unsigned char *) data;
+ tmp = cdat[0]; cdat[0] = cdat[1], cdat[1] = tmp;
+}*/
+
void Tqdcfr::FREADIA( unsigned num_ents, unsigned int* array ) {
unsigned rval = fread( array, sizeof(unsigned int), num_ents, cubFile );
IO_ASSERT( rval == num_ents );
+ if (swapForEndianness)
+ {
+ unsigned int * pt=array;
+ for (unsigned int i=0; i<num_ents; i++ )
+ {
+ swap4_uint((unsigned int *)pt);
+ pt++;
+ }
+ }
+
}
void Tqdcfr::FREADDA( unsigned num_ents, double* array ) {
unsigned rval = fread( array, sizeof(double), num_ents, cubFile );
IO_ASSERT( rval == num_ents );
+ if (swapForEndianness)
+ {
+ double * pt=array;
+ for (unsigned int i=0; i<num_ents; i++ )
+ {
+ swap8_voff((long *)pt);
+ pt++;
+ }
+ }
}
void Tqdcfr::FREADCA( unsigned num_ents, char* array ) {
@@ -172,7 +215,7 @@ ReaderIface* Tqdcfr::factory( Interface* iface )
Tqdcfr::Tqdcfr(Interface *impl)
: cubFile(NULL), globalIdTag(0), geomTag(0), uniqueIdTag(0),
blockTag(0), nsTag(0), ssTag(0), attribVectorTag(0), entityNameTag(0),
- categoryTag(0), hasMidNodesTag(0), printedSeqWarning(false), printedElemWarning(false)
+ categoryTag(0), hasMidNodesTag(0), swapForEndianness(false), printedSeqWarning(false), printedElemWarning(false)
{
assert(NULL != impl);
mdbImpl = impl;
@@ -1576,13 +1619,26 @@ ErrorCode Tqdcfr::read_file_header()
{
// read file header
FSEEK(4);
- FREADI(6);
- fileTOC.fileEndian = uint_buf[0];
- fileTOC.fileSchema = uint_buf[1];
- fileTOC.numModels = uint_buf[2];
- fileTOC.modelTableOffset = uint_buf[3];
- fileTOC.modelMetaDataOffset = uint_buf[4];
- fileTOC.activeFEModel = uint_buf[5];
+ // read tthe first int from the file
+ // if it is 0, it is littleEndian
+ unsigned rval = fread( &fileTOC.fileEndian, sizeof(unsigned int), 1, cubFile );
+ IO_ASSERT( rval == 1 );
+#ifdef WORDS_BIGENDIAN
+ if (fileTOC.fileEndian==0)
+ swapForEndianness=true;
+#else
+ if (fileTOC.fileEndian!=0)
+ swapForEndianness=true;
+#endif
+ if (debug)
+ std::cout << " swapping ? " << swapForEndianness << "\n";
+ FREADI(5);
+ //fileTOC.fileEndian = uint_buf[0];
+ fileTOC.fileSchema = uint_buf[0];
+ fileTOC.numModels = uint_buf[1];
+ fileTOC.modelTableOffset = uint_buf[2];
+ fileTOC.modelMetaDataOffset = uint_buf[3];
+ fileTOC.activeFEModel = uint_buf[4];
if (debug) fileTOC.print();
return MB_SUCCESS;
diff --git a/src/io/Tqdcfr.hpp b/src/io/Tqdcfr.hpp
index 1b4188f..eaf28ea 100644
--- a/src/io/Tqdcfr.hpp
+++ b/src/io/Tqdcfr.hpp
@@ -287,6 +287,7 @@ public:
attribVectorTag, entityNameTag, categoryTag, hasMidNodesTag;
std::map<int, EntityHandle> uidSetMap;
std::map<int, EntityHandle> gidSetMap[6];
+ bool swapForEndianness;
std::vector<unsigned int> uint_buf;
int *int_buf;
https://bitbucket.org/fathomteam/moab/commits/377b2a421b8f/
Changeset: 377b2a421b8f
Branch: None
User: danwu
Date: 2013-06-03 21:18:28
Summary: Update refactoring code for ReadNC class
Affected #: 11 files
diff --git a/src/io/Makefile.am b/src/io/Makefile.am
index 5a0a6f3..b7eae6a 100644
--- a/src/io/Makefile.am
+++ b/src/io/Makefile.am
@@ -18,6 +18,9 @@ if NETCDF_FILE
WriteSLAC.cpp WriteSLAC.hpp \
ReadNC.cpp ReadNC.hpp \
NCHelper.cpp NCHelper.hpp \
+ NCHelperEuler.cpp NCHelperEuler.hpp \
+ NCHelperFV.cpp NCHelperFV.hpp \
+ NCHelperHOMME.cpp NCHelperHOMME.hpp \
ReadGCRM.cpp ReadGCRM.hpp
else
MOAB_NETCDF_SRCS =
@@ -27,6 +30,9 @@ if PNETCDF_FILE
if !NETCDF_FILE
MOAB_NETCDF_SRCS += ReadNC.cpp ReadNC.hpp \
NCHelper.cpp NCHelper.hpp \
+ NCHelperEuler.cpp NCHelperEuler.hpp \
+ NCHelperFV.cpp NCHelperFV.hpp \
+ NCHelperHOMME.cpp NCHelperHOMME.hpp \
ReadGCRM.cpp ReadGCRM.hpp
endif
endif
diff --git a/src/io/NCHelper.cpp b/src/io/NCHelper.cpp
index d95cd34..cc38294 100644
--- a/src/io/NCHelper.cpp
+++ b/src/io/NCHelper.cpp
@@ -1,112 +1,45 @@
#include "NCHelper.hpp"
+#include "NCHelperEuler.hpp"
+#include "NCHelperFV.hpp"
+#include "NCHelperHOMME.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
namespace moab {
NCHelper* NCHelper::get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts)
{
- // Not a CAM file, will fill this in later for POP, CICE and CLM
- if (!readNC->isCam)
- return new (std::nothrow) NCHNotCam(readNC, fileId);
-
- // If a CAM file, which type?
- if (NCHEuler::can_read_file(readNC))
- return new (std::nothrow) NCHEuler(readNC, fileId);
- else if (NCHFV::can_read_file(readNC))
- return new (std::nothrow) NCHFV(readNC, fileId);
- else if (NCHHomme::can_read_file(readNC, opts))
- return new (std::nothrow) NCHHomme(readNC, fileId);
-
- // Unknown CAM grid
- return new (std::nothrow) NCHUnknownCam(readNC, fileId);
-}
-
-bool NCHEuler::can_read_file(ReadNC* readNC)
-{
- std::vector<std::string>& dimNames = readNC->dimNames;
-
- // If dimension names "lon" AND "lat' exist then it's the Eulerian Spectral grid or the FV grid
- if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
- dimNames.end(), std::string("lat")) != dimNames.end()))
- {
- // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
- if ((std::find(dimNames.begin(), dimNames.end(), std::string("slon")) != dimNames.end()) && (std::find(dimNames.begin(),
- dimNames.end(), std::string("slat")) != dimNames.end()))
- return false;
- else
- return true;
+ std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
+ if (attIt == readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
+ return NULL;
}
- return false;
-}
-
-ErrorCode NCHEuler::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
-{
- ErrorCode rval = _readNC->init_EulSpcscd_vals(opts, file_set);
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble initializing Euler grid.");
-
- return rval;
-}
-
-bool NCHFV::can_read_file(ReadNC* readNC)
-{
- std::vector<std::string>& dimNames = readNC->dimNames;
-
- // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
- if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
- dimNames.end(), std::string("lat")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slon"))
- != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end()))
- return true;
-
- return false;
-}
-
-ErrorCode NCHFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
-{
- ErrorCode rval = _readNC->init_FVCDscd_vals(opts, file_set);
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble initializing FV grid.");
-
- return rval;
-}
-
-bool NCHHomme::can_read_file(ReadNC* readNC, const FileOptions& opts)
-{
- // If global attribute "np" exists then it's the HOMME grid
- ErrorCode rval = readNC->check_np_attribute();
- if (MB_SUCCESS == rval) {
- if (MB_SUCCESS == opts.match_option("PARTITION_METHOD", "NODAL_PARTITION"))
- readNC->partMethod = -1;
-
- return true;
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ if (success != 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
+ return NULL;
}
- return false;
-}
-
-ErrorCode NCHHomme::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
-{
- ErrorCode rval = _readNC->init_HOMMEucd_vals();
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble initializing Homme grid.");
-
- return rval;
-}
-
-ErrorCode NCHUnknownCam::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
-{
- _readNC->readMeshIface->report_error("%s", "Unknown CAM grid.");
-
- return MB_FAILURE;
-}
-
-ErrorCode NCHNotCam::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
-{
- _readNC->readMeshIface->report_error("%s", "Unknown grid.");
+ // If a CAM file, which type?
+ if (att_data.find("CAM") != std::string::npos) {
+ int spectralOrder = -1;
+
+ if (NCHelperEuler::can_read_file(readNC))
+ return new (std::nothrow) NCHelperEuler(readNC, fileId);
+ else if (NCHelperFV::can_read_file(readNC))
+ return new (std::nothrow) NCHelperFV(readNC, fileId);
+ else if (NCHelperHOMME::can_read_file(readNC, opts, spectralOrder))
+ return new (std::nothrow) NCHelperHOMME(readNC, fileId, spectralOrder);
+ else // Unknown CAM type
+ return NULL;
+ }
- return MB_FAILURE;
+ // Not a CAM file, will fill this in later for POP, CICE and CLM
+ return NULL;
}
} // namespace moab
diff --git a/src/io/NCHelper.hpp b/src/io/NCHelper.hpp
index 10d6cdf..852fc61 100644
--- a/src/io/NCHelper.hpp
+++ b/src/io/NCHelper.hpp
@@ -22,137 +22,15 @@ public:
static NCHelper* get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts);
virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set) = 0;
-
virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads) = 0;
-
- virtual ErrorCode init_local_gid(EntityHandle file_set) = 0;
-
- virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums) = 0;
-
virtual std::string get_mesh_type_name() = 0;
+ virtual bool is_scd_mesh() = 0;
protected:
ReadNC* _readNC;
-
int _fileId;
};
-//! Child helper class for Eulerian Spectral grid (CAM_EUL)
-class NCHEuler : public NCHelper
-{
-public:
- NCHEuler(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
- static bool can_read_file(ReadNC* readNC);
-
-private:
- virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
-
- virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
- { return _readNC->create_scd_verts_quads(scdi, file_set, quads); }
-
- virtual ErrorCode init_local_gid(EntityHandle file_set)
- { return MB_SUCCESS; }
-
- virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
- { return _readNC->read_variables(file_set, var_names, tstep_nums, true); }
-
- virtual std::string get_mesh_type_name()
- { return "CAM_EUL"; }
-};
-
-//! Child helper class for Finite Volume grid (CAM_FV)
-class NCHFV : public NCHelper
-{
-public:
- NCHFV(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
- static bool can_read_file(ReadNC* readNC);
-
-private:
- virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
-
- virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
- { return _readNC->create_scd_verts_quads(scdi, file_set, quads); }
-
- virtual ErrorCode init_local_gid(EntityHandle file_set)
- { return MB_SUCCESS; }
-
- virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
- { return _readNC->read_variables(file_set, var_names, tstep_nums, true); }
-
- virtual std::string get_mesh_type_name()
- { return "CAM_FV"; }
-};
-
-//! Child helper class for HOMME grid (CAM_SE)
-class NCHHomme : public NCHelper
-{
-public:
- NCHHomme(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
- static bool can_read_file(ReadNC* readNC, const FileOptions& opts);
-
-private:
- virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
-
- virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
- { return _readNC->create_ucd_verts_quads(opts, file_set, quads); }
-
- virtual ErrorCode init_local_gid(EntityHandle file_set)
- { return _readNC->init_local_gid(file_set); }
-
- virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
- { return _readNC->read_variables(file_set, var_names, tstep_nums, false); }
-
- virtual std::string get_mesh_type_name()
- { return "CAM_SE"; }
-};
-
-//! Child helper class for unknown CAM grid (CAM_UNKNOWN)
-class NCHUnknownCam : public NCHelper
-{
-public:
- NCHUnknownCam(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
-private:
- virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
-
- virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
- { return MB_FAILURE; }
-
- virtual ErrorCode init_local_gid(EntityHandle file_set)
- { return MB_FAILURE; }
-
- virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
- { return MB_FAILURE; }
-
- virtual std::string get_mesh_type_name()
- { return "CAM_UNKNOWN"; }
-};
-
-//! Child helper class for unknown grid (NOT_CAM)
-class NCHNotCam : public NCHelper
-{
-public:
- NCHNotCam(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
-
-private:
- virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
-
- virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
- { return MB_FAILURE; }
-
- virtual ErrorCode init_local_gid(EntityHandle file_set)
- { return MB_FAILURE; }
-
- virtual ErrorCode read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
- { return MB_FAILURE; }
-
- virtual std::string get_mesh_type_name()
- { return "NOT_CAM"; }
-};
-
} // namespace moab
#endif
diff --git a/src/io/NCHelperEuler.cpp b/src/io/NCHelperEuler.cpp
new file mode 100644
index 0000000..2a6809d
--- /dev/null
+++ b/src/io/NCHelperEuler.cpp
@@ -0,0 +1,480 @@
+#include "NCHelperEuler.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+
+#include <cmath>
+#include <sstream>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) {_readNC->readMeshIface->report_error("%s", str); return rval;}
+
+#define ERRORS(err, str) \
+ if (err) {_readNC->readMeshIface->report_error("%s", str); return MB_FAILURE;}
+
+namespace moab {
+
+bool NCHelperEuler::can_read_file(ReadNC* readNC)
+{
+ std::vector<std::string>& dimNames = readNC->dimNames;
+
+ // If dimension names "lon" AND "lat' exist then it's the Eulerian Spectral grid or the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) != dimNames.end()))
+ {
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("slon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("slat")) != dimNames.end()))
+ return false;
+ else
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ ErrorCode rval = init_EulSpcscd_vals(opts, file_set);
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing Euler grid.");
+
+ return rval;
+}
+
+ErrorCode NCHelperEuler::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+{
+ ErrorCode rval = create_EulSpcscd_verts_quads(scdi, file_set, quads);
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble creating vertices and quads for Euler grid.");
+}
+
+ErrorCode NCHelperEuler::init_EulSpcscd_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ Interface*& mbImpl = _readNC->mbImpl;
+ std::vector<std::string>& dimNames = _readNC->dimNames;
+ std::vector<int>& dimVals = _readNC->dimVals;
+ std::string& iName = _readNC->iName;
+ std::string& jName = _readNC->jName;
+ std::string& tName = _readNC->tName;
+ std::string& iCName = _readNC->iCName;
+ std::string& jCName = _readNC->jCName;
+ std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
+ int& tMin = _readNC->tMin;
+ int& tMax = _readNC->tMax;
+ int (&gDims)[6] = _readNC->gDims;
+ int (&lDims)[6] = _readNC->lDims;
+ int (&lCDims)[6] = _readNC->lCDims;
+ int (&gCDims)[6] = _readNC->gCDims;
+ std::vector<double>& ilVals = _readNC->ilVals;
+ std::vector<double>& jlVals = _readNC->jlVals;
+ std::vector<double>& tVals = _readNC->tVals;
+ std::vector<double>& ilCVals = _readNC->ilCVals;
+ std::vector<double>& jlCVals = _readNC->jlCVals;
+ int& tDim = _readNC->tDim;
+ int& iCDim = _readNC->iCDim;
+ int& jCDim = _readNC->jCDim;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+ bool& isParallel = _readNC->isParallel;
+ int& partMethod = _readNC->partMethod;
+ int (&locallyPeriodic)[2] = _readNC->locallyPeriodic;
+ int (&globallyPeriodic)[2] = _readNC->globallyPeriodic;
+ ScdParData& parData = _readNC->parData;
+#ifdef USE_MPI
+ ParallelComm*& myPcomm = _readNC->myPcomm;
+#endif
+
+ // look for names of center i/j dimensions
+ std::vector<std::string>::iterator vit;
+ unsigned int idx;
+ iCName = std::string("lon");
+ iName = std::string("slon");
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), iCName.c_str())) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find center i variable.");
+ }
+ iCDim = idx;
+
+ // decide on i periodicity using math for now
+ std::vector<double> tilVals(dimVals[idx]);
+ ErrorCode rval = _readNC->read_coordinate(iCName.c_str(), 0, dimVals[idx] - 1, tilVals);
+ ERRORR(rval, "Trouble reading lon variable.");
+ if (std::fabs(2 * (*(tilVals.rbegin())) - *(tilVals.rbegin() + 1) - 360) < 0.001)
+ globallyPeriodic[0] = 1;
+
+ // now we can set gCDims and gDims for i
+ gCDims[0] = 0;
+ gDims[0] = 0;
+ gCDims[3] = dimVals[idx] - 1; // these are stored directly in file
+ gDims[3] = gCDims[3] + (globallyPeriodic[0] ? 0 : 1); // only if not periodic is vertex param max > elem param max
+
+ // now j
+ jCName = std::string("lat");
+ jName = std::string("slat");
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), jCName.c_str())) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find center j variable.");
+ }
+ jCDim = idx;
+
+ // for Eul models, will always be non-periodic in j
+ gCDims[1] = 0;
+ gDims[1] = 0;
+ gCDims[4] = dimVals[idx] - 1;
+ gDims[4] = gCDims[4] + 1;
+
+ // try a truly 2d mesh
+ gDims[2] = -1;
+ gDims[5] = -1;
+
+ // look for time dimensions
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ }
+ tDim = idx;
+ tMax = dimVals[idx] - 1;
+ tMin = 0;
+ tName = dimNames[idx];
+
+ // parse options to get subset
+ if (isParallel) {
+#ifdef USE_MPI
+ for (int i = 0; i < 6; i++)
+ parData.gDims[i] = gDims[i];
+ for (int i = 0; i < 3; i++)
+ parData.gPeriodic[i] = globallyPeriodic[i];
+ parData.partMethod = partMethod;
+ int pdims[3];
+
+ rval = ScdInterface::compute_partition(myPcomm->proc_config().proc_size(),
+ myPcomm->proc_config().proc_rank(),
+ parData, lDims, locallyPeriodic, pdims);
+ if (MB_SUCCESS != rval)
+ return rval;
+ for (int i = 0; i < 3; i++)
+ parData.pDims[i] = pdims[i];
+
+ dbgOut.tprintf(1, "Partition: %dx%d (out of %dx%d)\n",
+ lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
+ gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1);
+ if (myPcomm->proc_config().proc_rank() == 0)
+ dbgOut.tprintf(1, "Contiguous chunks of size %d bytes.\n", 8 * (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1));
+#endif
+ }
+ else {
+ for (int i = 0; i < 6; i++)
+ lDims[i] = gDims[i];
+ locallyPeriodic[0] = globallyPeriodic[0];
+ }
+
+ opts.get_int_option("IMIN", lDims[0]);
+ opts.get_int_option("IMAX", lDims[3]);
+ opts.get_int_option("JMIN", lDims[1]);
+ opts.get_int_option("JMAX", lDims[4]);
+
+ // now get actual coordinate values for vertices and cell centers; first resize
+ if (locallyPeriodic[0]) {
+ // if locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem coords
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0] + 1);
+ lCDims[3] = lDims[3];
+ }
+ else {
+ if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3]) {
+ // globally periodic and I'm the last proc, get fewer vertex coords than vertices in i
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] = lDims[3] - 1;
+ }
+ else {
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] = lDims[3] - 1;
+ }
+ }
+
+ lCDims[0] = lDims[0];
+ lCDims[4] = lDims[4] - 1;
+ lCDims[1] = lDims[1];
+
+ if (-1 != lDims[1]) {
+ jlVals.resize(lDims[4] - lDims[1] + 1);
+ jlCVals.resize(lCDims[4] - lCDims[1] + 1);
+ }
+
+ if (-1 != tMin)
+ tVals.resize(tMax - tMin + 1);
+
+ // now read coord values
+ std::map<std::string, ReadNC::VarData>::iterator vmit;
+ if (!ilCVals.empty()) {
+ if ((vmit = varInfo.find(iCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(iCName.c_str(), lDims[0], lDims[0] + ilCVals.size() - 1, ilCVals);
+ ERRORR(rval, "Trouble reading lon variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lon coordinate.");
+ }
+ }
+
+ if (!jlCVals.empty()) {
+ if ((vmit = varInfo.find(jCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(jCName.c_str(), lDims[1], lDims[1] + jlCVals.size() - 1, jlCVals);
+ ERRORR(rval, "Trouble reading lat variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lat coordinate.");
+ }
+ }
+
+ if (lDims[0] != -1) {
+ if ((vmit = varInfo.find(iCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ double dif = (ilCVals[1] - ilCVals[0]) / 2;
+ std::size_t i;
+ for (i = 0; i != ilCVals.size(); i++)
+ ilVals[i] = ilCVals[i] - dif;
+ // the last one is needed only if not periodic
+ if (!locallyPeriodic[0])
+ ilVals[i] = ilCVals[i - 1] + dif;
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ }
+ }
+
+ if (lDims[1] != -1) {
+ if ((vmit = varInfo.find(jCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ if (!isParallel || ((gDims[4] - gDims[1]) == (lDims[4] - lDims[1]))) {
+ std::string gwName("gw");
+ std::vector<double> gwVals(lDims[4] - lDims[1] - 1);
+ rval = _readNC->read_coordinate(gwName.c_str(), lDims[1], lDims[4] - 2, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ // copy the correct piece
+ jlVals[0] = -(M_PI / 2) * 180 / M_PI;
+ unsigned int i = 0;
+ double gwSum = -1;
+ for (i = 1; i != gwVals.size() + 1; i++) {
+ gwSum = gwSum + gwVals[i - 1];
+ jlVals[i] = std::asin(gwSum) * 180 / M_PI;
+ }
+ jlVals[i] = 90.0; // using value of i after loop exits.
+ }
+ else {
+ std::string gwName("gw");
+ double gwSum = 0;
+
+ // If this is the first row
+ if (lDims[1] == gDims[1]) {
+ std::vector<double> gwVals(lDims[4]);
+ rval = _readNC->read_coordinate(gwName.c_str(), 0, lDims[4] - 1, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ // copy the correct piece
+ jlVals[0] = -(M_PI / 2) * 180 / M_PI;
+ gwSum = -1;
+ for (std::size_t i = 1; i != jlVals.size(); i++) {
+ gwSum = gwSum + gwVals[i - 1];
+ jlVals[i] = std::asin(gwSum) * 180 / M_PI;
+ }
+ }
+ // or if it's the last row
+ else if (lDims[4] == gDims[4]) {
+ std::vector<double> gwVals(lDims[4] - 1);
+ rval = _readNC->read_coordinate(gwName.c_str(), 0, lDims[4] - 2, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ // copy the correct piece
+ gwSum = -1;
+ for (int j = 0; j != lDims[1] - 1; j++) {
+ gwSum = gwSum + gwVals[j];
+ }
+ std::size_t i = 0;
+ for (; i != jlVals.size() - 1; i++) {
+ gwSum = gwSum + gwVals[lDims[1] - 1 + i];
+ jlVals[i] = std::asin(gwSum) * 180 / M_PI;
+ }
+ jlVals[i] = 90.0; // using value of i after loop exits.
+ }
+ // it's in the middle
+ else {
+ int start = lDims[1] - 1;
+ int end = lDims[4] - 1;
+ std::vector<double> gwVals(end);
+ rval = _readNC->read_coordinate(gwName.c_str(), 0, end - 1, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ gwSum = -1;
+ for (int j = 0; j != start - 1; j++) {
+ gwSum = gwSum + gwVals[j];
+ }
+ std::size_t i = 0;
+ for (; i != jlVals.size(); i++) {
+ gwSum = gwSum + gwVals[start - 1 + i];
+ jlVals[i] = std::asin(gwSum) * 180 / M_PI;
+ }
+ }
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ }
+ }
+
+ if (tMin != -1) {
+ if ((vmit = varInfo.find(tName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+ }
+
+ dbgOut.tprintf(1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4]);
+ dbgOut.tprintf(1, "%d elements, %d vertices\n", (lDims[3] - lDims[0]) * (lDims[4] - lDims[1]), (lDims[3] - lDims[0] + 1)
+ * (lDims[4] - lDims[1] + 1));
+
+ // determine the entity location type of a variable
+ std::map<std::string, ReadNC::VarData>::iterator mit;
+ for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) {
+ ReadNC::VarData& vd = (*mit).second;
+ if ((std::find(vd.varDims.begin(), vd.varDims.end(), iCDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), jCDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCQUAD;
+ }
+
+ // <coordinate_dim_name>
+ std::vector<std::string> ijdimNames(4);
+ ijdimNames[0] = "__slon";
+ ijdimNames[1] = "__slat";
+ ijdimNames[2] = "__lon";
+ ijdimNames[3] = "__lat";
+
+ std::string tag_name;
+ int val_len = 0;
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ tag_name = ijdimNames[i];
+ void * val = NULL;
+ if (tag_name == "__slon") {
+ val = &ilVals[0];
+ val_len = ilVals.size();
+ }
+ else if (tag_name == "__slat") {
+ val = &jlVals[0];
+ val_len = jlVals.size();
+ }
+ else if (tag_name == "__lon") {
+ val = &ilCVals[0];
+ val_len = ilCVals.size();
+ }
+ else if (tag_name == "__lat") {
+ val = &jlCVals[0];
+ val_len = jlCVals.size();
+ }
+ Tag tagh = 0;
+ DataType data_type;
+
+ // assume all has same data type as lon
+ switch (varInfo["lon"].varDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ case NC_DOUBLE:
+ data_type = MB_TYPE_DOUBLE;
+ break;
+ case NC_FLOAT:
+ data_type = MB_TYPE_DOUBLE;
+ break;
+ case NC_INT:
+ data_type = MB_TYPE_INTEGER;
+ break;
+ case NC_SHORT:
+ default:
+ std::cerr << "Unrecognized data type for tag " << tag_name << std::endl;
+ ERRORR(MB_FAILURE, "Unrecognized data type");
+ break;
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, data_type, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating <coordinate_dim_name> tag.");
+ rval = mbImpl->tag_set_by_ptr(tagh, &file_set, 1, &val, &val_len);
+ ERRORR(rval, "Trouble setting data for <coordinate_dim_name> tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_LOC_MINMAX
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
+ tag_name = ss_tag_name.str();
+ Tag tagh = 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] == "__slon") {
+ val[0] = lDims[0];
+ val[1] = lDims[3];
+ }
+ else if (ijdimNames[i] == "__slat") {
+ val[0] = lDims[1];
+ val[1] = lDims[4];
+ }
+ else if (ijdimNames[i] == "__lon") {
+ val[0] = lCDims[0];
+ val[1] = lCDims[3];
+ }
+ else if (ijdimNames[i] == "__lat") {
+ val[0] = lCDims[1];
+ val[1] = lCDims[4];
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_LOC_MINMAX tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_LOC_MINMAX tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_GLOBAL_MINMAX
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
+ tag_name = ss_tag_name.str();
+ Tag tagh = 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] == "__slon") {
+ val[0] = gDims[0];
+ val[1] = gDims[3];
+ }
+ else if (ijdimNames[i] == "__slat") {
+ val[0] = gDims[1];
+ val[1] = gDims[4];
+ }
+ else if (ijdimNames[i] == "__lon") {
+ val[0] = gCDims[0];
+ val[1] = gCDims[3];
+ }
+ else if (ijdimNames[i] == "__lat") {
+ val[0] = gCDims[1];
+ val[1] = gCDims[4];
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_GLOBAL_MINMAX tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_GLOBAL_MINMAX tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // hack: create dummy tags, if needed, for variables like nbnd
+ // with no corresponding variables
+ _readNC->init_dims_with_no_cvars_info();
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelperEuler::create_EulSpcscd_verts_quads(ScdInterface* scdi, EntityHandle file_set, Range& quads)
+{
+ return _readNC->create_scd_verts_quads(scdi, file_set, quads);
+}
+
+} // namespace moab
diff --git a/src/io/NCHelperEuler.hpp b/src/io/NCHelperEuler.hpp
new file mode 100644
index 0000000..affea4a
--- /dev/null
+++ b/src/io/NCHelperEuler.hpp
@@ -0,0 +1,39 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelperEuler.hpp
+//
+// Purpose : Climate NC file helper for Eulerian Spectral grid
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPEREULER_HPP
+#define NCHELPEREULER_HPP
+
+#include "NCHelper.hpp"
+
+namespace moab {
+
+//! Child helper class for Eulerian Spectral grid (CAM_EUL)
+class NCHelperEuler : public NCHelper
+{
+public:
+ NCHelperEuler(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+ static bool can_read_file(ReadNC* readNC);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads);
+
+ virtual std::string get_mesh_type_name() { return "CAM_EUL"; }
+
+ virtual bool is_scd_mesh() { return true; }
+
+ ErrorCode init_EulSpcscd_vals(const FileOptions& opts, EntityHandle file_set);
+ ErrorCode create_EulSpcscd_verts_quads(ScdInterface* scdi, EntityHandle file_set, Range& quads);
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCHelperFV.cpp b/src/io/NCHelperFV.cpp
new file mode 100644
index 0000000..ff42fbe
--- /dev/null
+++ b/src/io/NCHelperFV.cpp
@@ -0,0 +1,476 @@
+#include "NCHelperFV.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+
+#include <cmath>
+#include <sstream>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) {_readNC->readMeshIface->report_error("%s", str); return rval;}
+
+namespace moab {
+
+bool NCHelperFV::can_read_file(ReadNC* readNC)
+{
+ std::vector<std::string>& dimNames = readNC->dimNames;
+
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slon"))
+ != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end()))
+ return true;
+
+ return false;
+}
+
+ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ ErrorCode rval = init_FVCDscd_vals(opts, file_set);
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing FV grid.");
+
+ return rval;
+}
+
+ErrorCode NCHelperFV::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+{
+ ErrorCode rval = create_FVCDscd_verts_quads(scdi, file_set, quads);
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble creating vertices and quads for FV grid.");
+
+ return rval;
+}
+
+ErrorCode NCHelperFV::init_FVCDscd_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ Interface*& mbImpl = _readNC->mbImpl;
+ std::vector<std::string>& dimNames = _readNC->dimNames;
+ std::vector<int>& dimVals = _readNC->dimVals;
+ std::string& iName = _readNC->iName;
+ std::string& jName = _readNC->jName;
+ std::string& tName = _readNC->tName;
+ std::string& iCName = _readNC->iCName;
+ std::string& jCName = _readNC->jCName;
+ std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
+ int& tMin = _readNC->tMin;
+ int& tMax = _readNC->tMax;
+ int (&gDims)[6] = _readNC->gDims;
+ int (&lDims)[6] = _readNC->lDims;
+ int (&gCDims)[6] = _readNC->gCDims;
+ int (&lCDims)[6] = _readNC->lCDims;
+ std::vector<double>& ilVals = _readNC->ilVals;
+ std::vector<double>& jlVals = _readNC->jlVals;
+ std::vector<double>& tVals = _readNC->tVals;
+ std::vector<double>& ilCVals = _readNC->ilCVals;
+ std::vector<double>& jlCVals = _readNC->jlCVals;
+ int& iDim = _readNC->iDim;
+ int& jDim = _readNC->jDim;
+ int& tDim = _readNC->tDim;
+ int& iCDim = _readNC->iCDim;
+ int& jCDim = _readNC->jCDim;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+ bool& isParallel = _readNC->isParallel;
+ int& partMethod = _readNC->partMethod;
+ int (&locallyPeriodic)[2] = _readNC->locallyPeriodic;
+ int (&globallyPeriodic)[2] = _readNC->globallyPeriodic;
+ ScdParData& parData = _readNC->parData;
+#ifdef USE_MPI
+ ParallelComm*& myPcomm = _readNC->myPcomm;
+#endif
+
+ std::vector<std::string>::iterator vit;
+ unsigned int idx;
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "slon")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find slon variable.");
+ }
+
+ iDim = idx;
+ gDims[3] = dimVals[idx] - 1;
+ gDims[0] = 0;
+ iName = dimNames[idx];
+
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "slat")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find slat variable.");
+ }
+ jDim = idx;
+ gDims[4] = dimVals[idx] - 1 + 2; // add 2 for the pole points
+ gDims[1] = 0;
+ jName = dimNames[idx];
+
+ // look for names of center i/j dimensions
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lon")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lon variable.");
+ }
+ iCDim = idx;
+ gCDims[3] = dimVals[idx] - 1;
+ gCDims[0] = 0;
+ iCName = dimNames[idx];
+
+ // check and set globallyPeriodic[0]
+ std::vector<double> til_vals(2);
+ ErrorCode rval = _readNC->read_coordinate(iCName.c_str(), dimVals[idx] - 2, dimVals[idx] - 1, til_vals);
+ ERRORR(rval, "Trouble reading slon variable.");
+ if (std::fabs(2 * til_vals[1] - til_vals[0] - 360) < 0.001)
+ globallyPeriodic[0] = 1;
+ if (globallyPeriodic[0])
+ assert("Number of vertices and edges should be same" && gDims[3] == gCDims[3]);
+ else
+ assert("Number of vertices should equal to number of edges plus one" && gDims[3] == gCDims[3]+1);
+
+#ifdef USE_MPI
+ // if serial, use a locally-periodic representation only if local mesh is periodic, otherwise don't
+ if ((isParallel && myPcomm->proc_config().proc_size() == 1) && globallyPeriodic[0])
+ locallyPeriodic[0] = 1;
+#else
+ if (globallyPeriodic[0])
+ locallyPeriodic[0] = 1;
+#endif
+
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lat")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lat variable.");
+ }
+ jCDim = idx;
+ gCDims[4] = dimVals[idx] - 1;
+ gCDims[1] = 0;
+ jCName = dimNames[idx];
+
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ }
+ tDim = idx;
+ tMax = dimVals[idx] - 1;
+ tMin = 0;
+ tName = dimNames[idx];
+
+ // parse options to get subset
+ if (isParallel) {
+#ifdef USE_MPI
+ for (int i = 0; i < 6; i++)
+ parData.gDims[i] = gDims[i];
+ for (int i = 0; i < 3; i++)
+ parData.gPeriodic[i] = globallyPeriodic[i];
+ parData.partMethod = partMethod;
+ int pdims[3];
+
+ rval = ScdInterface::compute_partition(myPcomm->proc_config().proc_size(),
+ myPcomm->proc_config().proc_rank(),
+ parData, lDims, locallyPeriodic, pdims);
+ if (MB_SUCCESS != rval)
+ return rval;
+ for (int i = 0; i < 3; i++)
+ parData.pDims[i] = pdims[i];
+
+ dbgOut.tprintf(1, "Partition: %dx%d (out of %dx%d)\n",
+ lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
+ gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1);
+ if (myPcomm->proc_config().proc_rank() == 0)
+ dbgOut.tprintf(1, "Contiguous chunks of size %d bytes.\n", 8 * (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1));
+#endif
+ }
+ else {
+ for (int i = 0; i < 6; i++)
+ lDims[i] = gDims[i];
+ locallyPeriodic[0] = globallyPeriodic[0];
+ }
+ opts.get_int_option("IMIN", lDims[0]);
+ opts.get_int_option("IMAX", lDims[3]);
+ opts.get_int_option("JMIN", lDims[1]);
+ opts.get_int_option("JMAX", lDims[4]);
+
+ // now get actual coordinate values for vertices and cell centers; first resize
+ if (locallyPeriodic[0]) {
+ // if locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem coords
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0] + 1);
+ lCDims[3] = lDims[3];
+ }
+ else {
+ if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3]) {
+ // globally periodic and I'm the last proc, get fewer vertex coords than vertices in i
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] = lDims[3] - 1;
+ }
+ else {
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] = lDims[3] - 1;
+ }
+ }
+
+ lCDims[0] = lDims[0];
+ lCDims[4] = lDims[4] - 1;
+ lCDims[1] = lDims[1];
+
+ if (-1 != lDims[1]) {
+ jlVals.resize(lDims[4] - lDims[1] + 1);
+ jlCVals.resize(lCDims[4] - lCDims[1] + 1);
+ }
+
+ if (-1 != tMin)
+ tVals.resize(tMax - tMin + 1);
+
+ // ... then read actual values
+ std::map<std::string, ReadNC::VarData>::iterator vmit;
+ if (lCDims[0] != -1) {
+ if ((vmit = varInfo.find(iCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(iCName.c_str(), lCDims[0], lCDims[3], ilCVals);
+ ERRORR(rval, "Trouble reading lon variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lon coordinate.");
+ }
+ }
+
+ if (lCDims[1] != -1) {
+ if ((vmit = varInfo.find(jCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(jCName.c_str(), lCDims[1], lCDims[4], jlCVals);
+ ERRORR(rval, "Trouble reading lat variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lat coordinate.");
+ }
+ }
+
+ if (lDims[0] != -1) {
+ if ((vmit = varInfo.find(iName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ // last column
+ if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3]) {
+ til_vals.resize(ilVals.size() - 1, 0.0);
+ rval = _readNC->read_coordinate(iName.c_str(), lDims[0], lDims[3] - 1, til_vals);
+ double dif = til_vals[1] - til_vals[0];
+ std::size_t i;
+ for (i = 0; i != til_vals.size(); i++)
+ ilVals[i] = til_vals[i];
+ ilVals[i] = ilVals[i - 1] + dif;
+ }
+ else {
+ rval = _readNC->read_coordinate(iName.c_str(), lDims[0], lDims[3], ilVals);
+ ERRORR(rval, "Trouble reading x variable.");
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ }
+ }
+
+ if (lDims[1] != -1) {
+ if ((vmit = varInfo.find(jName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ if (!isParallel || ((gDims[4] - gDims[1]) == (lDims[4] - lDims[1]))) {
+ std::vector<double> dummyVar(lDims[4] - lDims[1] - 1);
+ rval = _readNC->read_coordinate(jName.c_str(), lDims[1], lDims[4] - 2, dummyVar);
+ ERRORR(rval, "Trouble reading y variable.");
+ // copy the correct piece
+ jlVals[0] = -90.0;
+ unsigned int i = 0;
+ for (i = 1; i != dummyVar.size() + 1; i++)
+ jlVals[i] = dummyVar[i - 1];
+ jlVals[i] = 90.0; // using value of i after loop exits.
+ }
+ else {
+ // If this is the first row
+ // need to read one less then available and read it into a dummy var
+ if (lDims[1] == gDims[1]) {
+ std::vector<double> dummyVar(lDims[4] - lDims[1]);
+ rval = _readNC->read_coordinate(jName.c_str(), lDims[1], lDims[4] - 1, dummyVar);
+ ERRORR(rval, "Trouble reading y variable.");
+ // copy the correct piece
+ jlVals[0] = -90.0;
+ for (int i = 1; i < lDims[4] + 1; i++)
+ jlVals[i] = dummyVar[i - 1];
+ }
+ // or if it's the last row
+ else if (lDims[4] == gDims[4]) {
+ std::vector<double> dummyVar(lDims[4] - lDims[1]);
+ rval = _readNC->read_coordinate(jName.c_str(), lDims[1] - 1, lDims[4] - 2, dummyVar);
+ ERRORR(rval, "Trouble reading y variable.");
+ // copy the correct piece
+ std::size_t i = 0;
+ for (i = 0; i != dummyVar.size(); i++)
+ jlVals[i] = dummyVar[i];
+ jlVals[i] = 90.0; // using value of i after loop exits.
+ }
+ // it's in the middle
+ else {
+ rval = _readNC->read_coordinate(jCName.c_str(), lDims[1] - 1, lDims[4] - 1, jlVals);
+ ERRORR(rval, "Trouble reading y variable.");
+ }
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ }
+ }
+
+ if (tMin != -1) {
+ if ((vmit = varInfo.find(tName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+ }
+
+ dbgOut.tprintf(1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4]);
+ dbgOut.tprintf(1, "%d elements, %d vertices\n", (lDims[3] - lDims[0]) * (lDims[4] - lDims[1]), (lDims[3] - lDims[0] + 1)
+ * (lDims[4] - lDims[1] + 1));
+
+ // determine the entity location type of a variable
+ std::map<std::string, ReadNC::VarData>::iterator mit;
+ for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) {
+ ReadNC::VarData& vd = (*mit).second;
+ if ((std::find(vd.varDims.begin(), vd.varDims.end(), iCDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), jCDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCQUAD;
+ else if ((std::find(vd.varDims.begin(), vd.varDims.end(), jDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), iCDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCNSEDGE;
+ else if ((std::find(vd.varDims.begin(), vd.varDims.end(), jCDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), iDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCEWEDGE;
+ }
+
+ // <coordinate_dim_name>
+ std::vector<std::string> ijdimNames(4);
+ ijdimNames[0] = "__slon";
+ ijdimNames[1] = "__slat";
+ ijdimNames[2] = "__lon";
+ ijdimNames[3] = "__lat";
+
+ std::string tag_name;
+ int val_len = 0;
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ tag_name = ijdimNames[i];
+ void * val = NULL;
+ if (tag_name == "__slon") {
+ val = &ilVals[0];
+ val_len = ilVals.size();
+ }
+ else if (tag_name == "__slat") {
+ val = &jlVals[0];
+ val_len = jlVals.size();
+ }
+ else if (tag_name == "__lon") {
+ val = &ilCVals[0];
+ val_len = ilCVals.size();
+ }
+ else if (tag_name == "__lat") {
+ val = &jlCVals[0];
+ val_len = jlCVals.size();
+ }
+ Tag tagh = 0;
+ DataType data_type;
+
+ // assume all has same data type as lon
+ switch (varInfo["lon"].varDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ case NC_DOUBLE:
+ data_type = MB_TYPE_DOUBLE;
+ break;
+ case NC_FLOAT:
+ data_type = MB_TYPE_DOUBLE;
+ break;
+ case NC_INT:
+ data_type = MB_TYPE_INTEGER;
+ break;
+ case NC_SHORT:
+ default:
+ std::cerr << "Unrecognized data type for tag " << tag_name << std::endl;
+ ERRORR(MB_FAILURE, "Unrecognized data type");
+ break;
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, data_type, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating <coordinate_dim_name> tag.");
+ rval = mbImpl->tag_set_by_ptr(tagh, &file_set, 1, &val, &val_len);
+ ERRORR(rval, "Trouble setting data for <coordinate_dim_name> tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_LOC_MINMAX
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
+ tag_name = ss_tag_name.str();
+ Tag tagh = 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] == "__slon") {
+ val[0] = lDims[0];
+ val[1] = lDims[3];
+ }
+ else if (ijdimNames[i] == "__slat") {
+ val[0] = lDims[1];
+ val[1] = lDims[4];
+ }
+ else if (ijdimNames[i] == "__lon") {
+ val[0] = lCDims[0];
+ val[1] = lCDims[3];
+ }
+ else if (ijdimNames[i] == "__lat") {
+ val[0] = lCDims[1];
+ val[1] = lCDims[4];
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_LOC_MINMAX tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_LOC_MINMAX tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_GLOBAL_MINMAX
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
+ tag_name = ss_tag_name.str();
+ Tag tagh = 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] == "__slon") {
+ val[0] = gDims[0];
+ val[1] = gDims[3];
+ }
+ else if (ijdimNames[i] == "__slat") {
+ val[0] = gDims[1];
+ val[1] = gDims[4];
+ }
+ else if (ijdimNames[i] == "__lon") {
+ val[0] = gCDims[0];
+ val[1] = gCDims[3];
+ }
+ else if (ijdimNames[i] == "__lat") {
+ val[0] = gCDims[1];
+ val[1] = gCDims[4];
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_GLOBAL_MINMAX tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_GLOBAL_MINMAX tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // hack: create dummy tags, if needed, for dimensions like nbnd
+ // with no corresponding variables
+ _readNC->init_dims_with_no_cvars_info();
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelperFV::create_FVCDscd_verts_quads(ScdInterface* scdi, EntityHandle file_set, Range& quads)
+{
+ return _readNC->create_scd_verts_quads(scdi, file_set, quads);
+}
+
+} // namespace moab
diff --git a/src/io/NCHelperFV.hpp b/src/io/NCHelperFV.hpp
new file mode 100644
index 0000000..b9f5d31
--- /dev/null
+++ b/src/io/NCHelperFV.hpp
@@ -0,0 +1,36 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelperFV.hpp
+//
+// Purpose : Climate NC file helper for Finite Volume grid
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPERFV_HPP
+#define NCHELPERFV_HPP
+
+#include "NCHelper.hpp"
+
+namespace moab {
+
+//! Child helper class for Finite Volume grid (CAM_FV)
+class NCHelperFV : public NCHelper
+{
+public:
+ NCHelperFV(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+ static bool can_read_file(ReadNC* readNC);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads);
+ virtual std::string get_mesh_type_name() { return "CAM_FV"; }
+ virtual bool is_scd_mesh() { return true; }
+
+private:
+ ErrorCode init_FVCDscd_vals(const FileOptions& opts, EntityHandle file_set);
+ ErrorCode create_FVCDscd_verts_quads(ScdInterface* scdi, EntityHandle file_set, Range& quads);
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCHelperHOMME.cpp b/src/io/NCHelperHOMME.cpp
new file mode 100644
index 0000000..604e000
--- /dev/null
+++ b/src/io/NCHelperHOMME.cpp
@@ -0,0 +1,501 @@
+#include "NCHelperHOMME.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+#include "moab/SpectralMeshTool.hpp"
+
+#include <cmath>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) {_readNC->readMeshIface->report_error("%s", str); return rval;}
+
+#define ERRORS(err, str) \
+ if (err) {_readNC->readMeshIface->report_error("%s", str); return MB_FAILURE;}
+
+namespace moab {
+
+bool NCHelperHOMME::can_read_file(ReadNC* readNC, const FileOptions& opts, int& spectralOrder)
+{
+ // If global attribute "np" exists then it's the HOMME grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("np");
+ if (attIt != readNC->globalAtts.end())
+ {
+ int success = NCFUNC(get_att_int)(readNC->fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &spectralOrder);
+ if (success != 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read np global attribute int data.");
+ return false;
+ }
+
+ spectralOrder--; // Spectral order is one less than np
+
+ if (MB_SUCCESS == opts.match_option("PARTITION_METHOD", "NODAL_PARTITION"))
+ readNC->partMethod = -1;
+
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHelperHOMME::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ ErrorCode rval = init_HOMMEucd_vals();
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble initializing HOMME grid.");
+
+ return rval;
+}
+
+ErrorCode NCHelperHOMME::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+{
+ ErrorCode rval = create_HOMMEucd_verts_quads(opts, file_set, quads);
+ if (MB_SUCCESS != rval)
+ _readNC->readMeshIface->report_error("%s", "Trouble creating vertices and quads for HOMME grid.");
+
+ return rval;
+}
+
+ErrorCode NCHelperHOMME::init_HOMMEucd_vals()
+{
+ std::vector<std::string>& dimNames = _readNC->dimNames;
+ std::vector<int>& dimVals = _readNC->dimVals;
+ std::string& kName = _readNC->kName;
+ std::string& tName = _readNC->tName;
+ std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
+ int& tMin = _readNC->tMin;
+ int& tMax = _readNC->tMax;
+ int (&gDims)[6] = _readNC->gDims;
+ int (&lDims)[6] = _readNC->lDims;
+ int& iDim = _readNC->iDim;
+ int& kDim = _readNC->kDim;
+ int& tDim = _readNC->tDim;
+ std::vector<double>& ilVals = _readNC->ilVals;
+ std::vector<double>& jlVals = _readNC->jlVals;
+ std::vector<double>& klVals = _readNC->klVals;
+ std::vector<double>& tVals = _readNC->tVals;
+
+ ErrorCode rval;
+ unsigned int idx;
+ std::vector<std::string>::iterator vit;
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ tDim = idx;
+ tMax = dimVals[idx] - 1;
+ tMin = 0;
+ tName = dimNames[idx];
+
+ // get number of vertices (labeled as number of columns) and levels
+ gDims[0] = gDims[3] = -1;
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "ncol")) != dimNames.end()) {
+ idx = vit - dimNames.begin();
+ gDims[3] = dimVals[idx] - 1;
+ gDims[0] = 0;
+ iDim = idx;
+ }
+ if (-1 == gDims[0])
+ return MB_FAILURE;
+
+ // set j coordinate to the number of quads
+ gDims[1] = gDims[0];
+ gDims[4] = gDims[3] - 2;
+
+ gDims[2] = gDims[5] = -1;
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end()) {
+ idx = vit - dimNames.begin();
+ gDims[5] = dimVals[idx] - 1, gDims[2] = 0, kName = std::string("lev");
+ kDim = idx;
+ }
+ if (-1 == gDims[2])
+ return MB_FAILURE;
+
+ // read coordinate data
+ std::map<std::string, ReadNC::VarData>::iterator vmit;
+ if (gDims[0] != -1) {
+ if ((vmit = varInfo.find("lon")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate("lon", gDims[0], gDims[3], ilVals);
+ ERRORR(rval, "Trouble reading x variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ }
+ }
+
+ // store lat values in jlVals parameterized by j
+ if (gDims[1] != -1) {
+ if ((vmit = varInfo.find("lat")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate("lat", gDims[0], gDims[3], jlVals);
+ ERRORR(rval, "Trouble reading y variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ }
+ }
+
+ if (gDims[2] != -1) {
+ if ((vmit = varInfo.find("lev")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate("lev", gDims[2], gDims[5], klVals);
+ ERRORR(rval, "Trouble reading z variable.");
+
+ // decide whether down is positive
+ char posval[10];
+ int success = NCFUNC(get_att_text)(_readNC->fileId, (*vmit).second.varId, "positive", posval);
+ if (0 == success && !strcmp(posval, "down")) {
+ for (std::vector<double>::iterator dvit = klVals.begin(); dvit != klVals.end(); ++dvit)
+ (*dvit) *= -1.0;
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find z coordinate.");
+ }
+ }
+
+ if (tMin != -1) {
+ if ((vmit = varInfo.find(tName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+ }
+
+ if ((vmit = varInfo.find(tName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+
+ // determine the entity location type of a variable
+ std::map<std::string, ReadNC::VarData>::iterator mit;
+ for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) {
+ ReadNC::VarData& vd = (*mit).second;
+ if ((std::find(vd.varDims.begin(), vd.varDims.end(), iDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), kDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCNODE;
+ }
+
+ std::copy(gDims, gDims + 6, lDims);
+
+ // don't read coordinates of columns until we actually create the mesh
+
+ // hack: create dummy tags, if needed, for variables like ncol and nbnd
+ // with no corresponding variables
+ _readNC->init_dims_with_no_cvars_info();
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelperHOMME::create_HOMMEucd_verts_quads(const FileOptions& opts, EntityHandle file_set, Range& quads)
+{
+ Interface*& mbImpl = _readNC->mbImpl;
+ std::string& fileName = _readNC->fileName;
+ int& connectId = _readNC->connectId;
+ int (&gDims)[6] = _readNC->gDims;
+ int (&lDims)[6] = _readNC->lDims;
+ std::vector<double>& ilVals = _readNC->ilVals;
+ std::vector<double>& jlVals = _readNC->jlVals;
+ std::vector<double>& klVals = _readNC->klVals;
+ Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
+ const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+ bool& isParallel = _readNC->isParallel;
+ Range& localGid = _readNC->localGid;
+#ifdef USE_MPI
+ ParallelComm*& myPcomm = _readNC->myPcomm;
+#endif
+ bool& spectralMesh = _readNC->spectralMesh;
+
+ // need to get/read connectivity data before creating elements
+ std::string conn_fname;
+
+ // try to open the connectivity file through CONN option, if used
+ ErrorCode rval = opts.get_str_option("CONN", conn_fname);
+ if (MB_SUCCESS != rval) {
+ // default convention for reading HOMME is a file HommeMapping.nc in same dir as data file
+ conn_fname = std::string(fileName);
+ size_t idx = conn_fname.find_last_of("/");
+ if (idx != std::string::npos)
+ conn_fname = conn_fname.substr(0, idx).append("/HommeMapping.nc");
+ else
+ conn_fname = "HommeMapping.nc";
+ }
+
+ int success;
+
+ int rank, procs;
+#ifdef PNETCDF_FILE
+ if (isParallel) {
+ success = NCFUNC(open)(myPcomm->proc_config().proc_comm(), conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId);
+ rank = myPcomm->proc_config().proc_rank();
+ procs = myPcomm->proc_config().proc_size();
+ }
+ else {
+ success = NCFUNC(open)(MPI_COMM_SELF, conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId);
+ rank = 0;
+ procs = 1;
+ }
+#else
+ success = NCFUNC(open)(conn_fname.c_str(), 0, &connectId);
+ rank = 0;
+ procs = 1;
+#endif
+ ERRORS(success, "Failed on open.");
+
+ std::vector<std::string> conn_names;
+ std::vector<int> conn_vals;
+ rval = _readNC->get_dimensions(connectId, conn_names, conn_vals);
+ ERRORR(rval, "Failed to get dimensions for connectivity.");
+
+ if (conn_vals[0] != gDims[3] - gDims[0] + 1 - 2) {
+ dbgOut.tprintf(1, "Warning: number of quads from %s and vertices from %s are inconsistent; nverts = %d, nquads = %d.\n",
+ conn_fname.c_str(), fileName.c_str(), gDims[3] - gDims[0] + 1, conn_vals[0]);
+ }
+
+ // read connectivity into temporary variable
+ int num_fine_quads, num_coarse_quads, start_idx;
+ std::vector<std::string>::iterator vit;
+ int idx;
+ if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncells")) != conn_names.end())
+ {
+ idx = vit - conn_names.begin();
+ }
+ else if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncenters")) != conn_names.end())
+ {
+ idx = vit - conn_names.begin();
+ }
+ else
+ {
+ ERRORR(MB_FAILURE, "Failed to get number of quads.");
+ }
+ int num_quads = conn_vals[idx];
+
+ // get the connectivity into tmp_conn2 and permute into tmp_conn
+ int cornerVarId;
+ success = NCFUNC(inq_varid)(connectId, "element_corners", &cornerVarId);
+ ERRORS(success, "Failed to get variable id.");
+ NCDF_SIZE tmp_dims[2] = { 0, 0 }, tmp_counts[2] = { 4, static_cast<size_t>(num_quads) };
+ std::vector<int> tmp_conn(4 * num_quads), tmp_conn2(4 * num_quads);
+ success = NCFUNCAG(_vara_int)(connectId, cornerVarId, tmp_dims, tmp_counts, &tmp_conn2[0] NCREQ);
+ ERRORS(success, "Failed to get temporary connectivity.");
+ success = NCFUNC(close)(connectId);
+ ERRORS(success, "Failed on close.");
+ // permute the connectivity
+ for (int i = 0; i < num_quads; i++) {
+ tmp_conn[4*i] = tmp_conn2[i];
+ tmp_conn[4*i + 1] = tmp_conn2[i + 1 * num_quads];
+ tmp_conn[4*i + 2] = tmp_conn2[i + 2 * num_quads];
+ tmp_conn[4*i + 3] = tmp_conn2[i + 3 * num_quads];
+ }
+
+ // need to know whether we'll be creating gather mesh later, to make sure we allocate enough space
+ // in one shot
+ bool create_gathers = true;
+#ifdef USE_MPI
+ if (isParallel)
+ if (myPcomm->proc_config().proc_rank() != 0)
+ create_gathers = false;
+#endif
+
+ // compute the number of local quads, accounting for coarse or fine representation
+ // spectral_unit is the # fine quads per coarse quad, or spectralOrder^2
+ int spectral_unit = (spectralMesh ? _spectralOrder * _spectralOrder : 1);
+ // num_coarse_quads is the number of quads instantiated in MOAB; if !spectralMesh, num_coarse_quads = num_fine_quads
+ num_coarse_quads = int(std::floor(1.0 * num_quads / (spectral_unit * procs)));
+ // start_idx is the starting index in the HommeMapping connectivity list for this proc, before converting to coarse quad representation
+ start_idx = 4 * rank * num_coarse_quads * spectral_unit;
+ // iextra = # coarse quads extra after equal split over procs
+ int iextra = num_quads % (procs * spectral_unit);
+ if (rank < iextra) num_coarse_quads++;
+ start_idx += 4 * spectral_unit * std::min(rank, iextra);
+ // num_fine_quads is the number of quads in the connectivity list in HommeMapping file assigned to this proc
+ num_fine_quads = spectral_unit * num_coarse_quads;
+
+ // now create num_coarse_quads
+ EntityHandle *conn_arr;
+ EntityHandle start_vertex;
+ Range tmp_range;
+
+ // read connectivity into that space
+ EntityHandle *sv_ptr = NULL, start_quad;
+ SpectralMeshTool smt(mbImpl, _spectralOrder);
+ if (!spectralMesh) {
+ rval = _readNC->readMeshIface->get_element_connect(num_coarse_quads, 4,
+ MBQUAD, 0, start_quad, conn_arr,
+ // might have to create gather mesh later
+ (create_gathers ? num_coarse_quads + num_quads : num_coarse_quads));
+ ERRORR(rval, "Failed to create quads.");
+ tmp_range.insert(start_quad, start_quad + num_coarse_quads - 1);
+ std::copy(&tmp_conn[start_idx], &tmp_conn[start_idx + 4 * num_fine_quads], conn_arr);
+ std::copy(conn_arr, conn_arr + 4 * num_fine_quads, range_inserter(localGid));
+ }
+ else {
+ rval = smt.create_spectral_elems(&tmp_conn[0], num_fine_quads, 2, tmp_range, start_idx, &localGid);
+ ERRORR(rval, "Failed to create spectral elements.");
+ int count, v_per_e;
+ rval = mbImpl->connect_iterate(tmp_range.begin(), tmp_range.end(), conn_arr, v_per_e, count);
+ ERRORR(rval, "Failed to get connectivity of spectral elements.");
+ rval = mbImpl->tag_iterate(smt.spectral_vertices_tag(true), tmp_range.begin(), tmp_range.end(),
+ count, (void*&)sv_ptr);
+ ERRORR(rval, "Failed to get fine connectivity of spectral elements.");
+ }
+
+ // on this proc, I get columns ldims[1]..ldims[4], inclusive; need to find which vertices those correpond to
+ unsigned int num_local_verts = localGid.size();
+ unsigned int num_total_verts = gDims[3] - gDims[0] + 1;
+
+ // create vertices
+ std::vector<double*> arrays;
+ rval = _readNC->readMeshIface->get_node_coords(3, num_local_verts, 0, start_vertex, arrays,
+ // might have to create gather mesh later
+ (create_gathers ? num_local_verts+num_total_verts : num_local_verts));
+ ERRORR(rval, "Couldn't create vertices in ucd mesh.");
+
+ // set vertex coordinates
+ Range::iterator rit;
+ double *xptr = arrays[0], *yptr = arrays[1], *zptr = arrays[2];
+ int i;
+ for (i = 0, rit = localGid.begin(); i < (int)num_local_verts; i++, ++rit) {
+ assert(*rit < ilVals.size()+1);
+ xptr[i] = ilVals[(*rit) - 1];
+ yptr[i] = jlVals[(*rit) - 1];
+ zptr[i] = klVals[lDims[2]];
+ }
+
+ const double pideg = acos(-1.0) / 180.0;
+ for (i = 0; i < (int)num_local_verts; i++) {
+ double cosphi = cos(pideg * yptr[i]);
+ double zmult = sin(pideg * yptr[i]), xmult = cosphi * cos(xptr[i] * pideg), ymult = cosphi * sin(xptr[i] * pideg);
+ double rad = 8.0e3 + klVals[lDims[2]];
+ xptr[i] = rad * xmult, yptr[i] = rad * ymult, zptr[i] = rad * zmult;
+ }
+
+ // get ptr to gid memory for vertices
+ Range vert_range(start_vertex, start_vertex + num_local_verts - 1);
+ void *data;
+ int count;
+ rval = mbImpl->tag_iterate(mGlobalIdTag, vert_range.begin(), vert_range.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator.");
+ assert(count == (int) num_local_verts);
+ int *gid_data = (int*) data;
+ std::copy(localGid.begin(), localGid.end(), gid_data);
+ // duplicate global id data, which will be used to resolve sharing
+ if (mpFileIdTag)
+ {
+ rval = mbImpl->tag_iterate(*mpFileIdTag, vert_range.begin(), vert_range.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator on file id tag.");
+ assert(count == (int) num_local_verts);
+ gid_data = (int*) data;
+ std::copy(localGid.begin(), localGid.end(), gid_data);
+ }
+
+ // create map from file ids to vertex handles, used later to set connectivity
+ std::map<EntityHandle, EntityHandle> vert_handles;
+ for (rit = localGid.begin(), i = 0; rit != localGid.end(); ++rit, i++) {
+ vert_handles[*rit] = start_vertex + i;
+ }
+
+ // compute proper handles in connectivity using offset
+ for (int q = 0; q < 4 * num_coarse_quads; q++) {
+ conn_arr[q] = vert_handles[conn_arr[q]];
+ assert(conn_arr[q]);
+ }
+ if (spectralMesh) {
+ int verts_per_quad = (_spectralOrder + 1) * (_spectralOrder + 1);
+ for (int q = 0; q < verts_per_quad * num_coarse_quads; q++) {
+ sv_ptr[q] = vert_handles[sv_ptr[q]];
+ assert(sv_ptr[q]);
+ }
+ }
+
+ // add new vertices and elements to the set
+ quads.merge(tmp_range);
+ tmp_range.insert(start_vertex, start_vertex + num_local_verts - 1);
+ rval = mbImpl->add_entities(file_set, tmp_range);
+ ERRORR(rval, "Couldn't add new vertices and quads/hexes to file set.");
+
+ // mark the set with the spectral order
+ Tag sporder;
+ rval = mbImpl->tag_get_handle("SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, sporder, MB_TAG_CREAT | MB_TAG_SPARSE);
+ ERRORR(rval, "Couldn't create spectral order tag.");
+ rval = mbImpl->tag_set_data(sporder, &file_set, 1, &_spectralOrder);
+ ERRORR(rval, "Couldn't set value for spectral order tag.");
+
+#ifdef USE_MPI
+ if (isParallel && myPcomm->proc_config().proc_rank() == 0) {
+#endif
+ EntityHandle gather_set;
+ rval = mbImpl->create_meshset(MESHSET_SET, gather_set);
+ ERRORR(rval, "Trouble creating gather set.");
+
+ // create vertices
+ arrays.clear();
+ // don't need to specify allocation number here, because we know enough verts were created before
+ rval = _readNC->readMeshIface->get_node_coords(3, num_total_verts, 0, start_vertex, arrays);
+ ERRORR(rval, "Couldn't create vertices in ucd mesh for gather set.");
+
+ xptr = arrays[0], yptr = arrays[1], zptr = arrays[2];
+ for (i = 0; i < (int)num_total_verts; i++) {
+ double cosphi = cos(pideg * jlVals[i]);
+ double zmult = sin(pideg * jlVals[i]);
+ double xmult = cosphi * cos(ilVals[i] * pideg);
+ double ymult = cosphi * sin(ilVals[i] * pideg);
+ double rad = 8.0e3 + klVals[lDims[2]];
+ xptr[i] = rad * xmult;
+ yptr[i] = rad * ymult;
+ zptr[i] = rad * zmult;
+ }
+
+ // get ptr to gid memory for vertices
+ Range gather_verts(start_vertex, start_vertex + num_total_verts - 1);
+ rval = mbImpl->tag_iterate(mGlobalIdTag, gather_verts.begin(), gather_verts.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator.");
+ assert(count == (int) num_total_verts);
+ gid_data = (int*) data;
+ for (int j = 1; j <= (int) num_total_verts; j++)
+ gid_data[j - 1] = j;
+ // set the file id tag too, it should be bigger something not interfering with global id
+ if (mpFileIdTag) {
+ rval = mbImpl->tag_iterate(*mpFileIdTag, gather_verts.begin(), gather_verts.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator in file id tag.");
+ assert(count == (int) num_total_verts);
+ gid_data = (int*) data;
+ for (int j = 1; j <= (int) num_total_verts; j++)
+ gid_data[j - 1] = num_total_verts + j; // bigger than global id tag
+ }
+
+ rval = mbImpl->add_entities(gather_set, gather_verts);
+ ERRORR(rval, "Couldn't add vertices to gather set.");
+
+ // create quads
+ Range gather_quads;
+ // don't need to specify allocation number here, because we know enough quads were created before
+ rval = _readNC->readMeshIface->get_element_connect(num_quads, 4,
+ MBQUAD, 0, start_quad, conn_arr);
+ ERRORR(rval, "Failed to create quads.");
+ gather_quads.insert(start_quad, start_quad + num_quads - 1);
+ std::copy(&tmp_conn[0], &tmp_conn[4 * num_quads], conn_arr);
+ for (i = 0; i != 4 * num_quads; i++)
+ conn_arr[i] += start_vertex - 1; // connectivity array is shifted by where the gather verts start
+ rval = mbImpl->add_entities(gather_set, gather_quads);
+ ERRORR(rval, "Couldn't add quads to gather set.");
+
+ Tag gathersettag;
+ rval = mbImpl->tag_get_handle("GATHER_SET", 1, MB_TYPE_INTEGER, gathersettag,
+ MB_TAG_CREAT | MB_TAG_SPARSE);
+ ERRORR(rval, "Couldn't create gather set tag.");
+ int gatherval = 1;
+ rval = mbImpl->tag_set_data(gathersettag, &gather_set, 1, &gatherval);
+ ERRORR(rval, "Couldn't set value for gather set tag.");
+
+#ifdef USE_MPI
+ }
+#endif
+
+ return MB_SUCCESS;
+}
+
+} // namespace moab
diff --git a/src/io/NCHelperHOMME.hpp b/src/io/NCHelperHOMME.hpp
new file mode 100644
index 0000000..f0dc1cd
--- /dev/null
+++ b/src/io/NCHelperHOMME.hpp
@@ -0,0 +1,38 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelperHOMME.hpp
+//
+// Purpose : Climate NC file helper for HOMME grid
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPERHOMME_HPP
+#define NCHELPERHOMME_HPP
+
+#include "NCHelper.hpp"
+
+namespace moab {
+
+//! Child helper class for HOMME grid (CAM_SE)
+class NCHelperHOMME : public NCHelper
+{
+public:
+ NCHelperHOMME(ReadNC* readNC, int fileId, int spectralOrder) : NCHelper(readNC, fileId), _spectralOrder(spectralOrder) {}
+ static bool can_read_file(ReadNC* readNC, const FileOptions& opts, int& spectralOrder);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads);
+ virtual std::string get_mesh_type_name() { return "CAM_SE"; }
+ virtual bool is_scd_mesh() { return false; }
+
+ ErrorCode init_HOMMEucd_vals();
+ ErrorCode create_HOMMEucd_verts_quads(const FileOptions& opts, EntityHandle file_set, Range& quads);
+
+private:
+ int _spectralOrder; // read from variable 'np'
+};
+
+} // namespace moab
+
+#endif
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/fathomteam/moab/commits/9fa973ff0484/
Changeset: 9fa973ff0484
Branch: None
User: danwu
Date: 2013-06-05 16:56:37
Summary: Update refactoring code for ReadNC class
Affected #: 9 files
diff --git a/src/io/NCHelper.cpp b/src/io/NCHelper.cpp
index cc38294..e0f021e 100644
--- a/src/io/NCHelper.cpp
+++ b/src/io/NCHelper.cpp
@@ -8,38 +8,14 @@ namespace moab {
NCHelper* NCHelper::get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts)
{
- std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
- if (attIt == readNC->globalAtts.end()) {
- readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
+ if (NCHelperEuler::can_read_file(readNC, fileId))
+ return new (std::nothrow) NCHelperEuler(readNC, fileId);
+ else if (NCHelperFV::can_read_file(readNC, fileId))
+ return new (std::nothrow) NCHelperFV(readNC, fileId);
+ else if (NCHelperHOMME::can_read_file(readNC, fileId))
+ return new (std::nothrow) NCHelperHOMME(readNC, fileId, opts);
+ else // Unknown NetCDF grid (will fill this in later for POP, CICE and CLM)
return NULL;
- }
-
- unsigned int sz = attIt->second.attLen;
- std::string att_data;
- att_data.resize(sz + 1);
- att_data[sz] = '\000';
- int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
- if (success != 0) {
- readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
- return NULL;
- }
-
- // If a CAM file, which type?
- if (att_data.find("CAM") != std::string::npos) {
- int spectralOrder = -1;
-
- if (NCHelperEuler::can_read_file(readNC))
- return new (std::nothrow) NCHelperEuler(readNC, fileId);
- else if (NCHelperFV::can_read_file(readNC))
- return new (std::nothrow) NCHelperFV(readNC, fileId);
- else if (NCHelperHOMME::can_read_file(readNC, opts, spectralOrder))
- return new (std::nothrow) NCHelperHOMME(readNC, fileId, spectralOrder);
- else // Unknown CAM type
- return NULL;
- }
-
- // Not a CAM file, will fill this in later for POP, CICE and CLM
- return NULL;
}
} // namespace moab
diff --git a/src/io/NCHelperEuler.cpp b/src/io/NCHelperEuler.cpp
index 2a6809d..9346d22 100644
--- a/src/io/NCHelperEuler.cpp
+++ b/src/io/NCHelperEuler.cpp
@@ -13,20 +13,37 @@
namespace moab {
-bool NCHelperEuler::can_read_file(ReadNC* readNC)
+bool NCHelperEuler::can_read_file(ReadNC* readNC, int fileId)
{
std::vector<std::string>& dimNames = readNC->dimNames;
- // If dimension names "lon" AND "lat' exist then it's the Eulerian Spectral grid or the FV grid
+ // If dimension names "lon" AND "lat' exist then it could be either the Eulerian Spectral grid or the FV grid
if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
- dimNames.end(), std::string("lat")) != dimNames.end()))
- {
- // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
+ dimNames.end(), std::string("lat")) != dimNames.end())) {
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it should be the FV grid
if ((std::find(dimNames.begin(), dimNames.end(), std::string("slon")) != dimNames.end()) && (std::find(dimNames.begin(),
dimNames.end(), std::string("slat")) != dimNames.end()))
return false;
- else
- return true;
+
+ // Make sure it is CAM grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
+ if (attIt == readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
+ return false;
+ }
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ if (success != 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
+ return false;
+ }
+ if (att_data.find("CAM") == std::string::npos)
+ return false;
+
+ return true;
}
return false;
@@ -34,22 +51,6 @@ bool NCHelperEuler::can_read_file(ReadNC* readNC)
ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
{
- ErrorCode rval = init_EulSpcscd_vals(opts, file_set);
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble initializing Euler grid.");
-
- return rval;
-}
-
-ErrorCode NCHelperEuler::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
-{
- ErrorCode rval = create_EulSpcscd_verts_quads(scdi, file_set, quads);
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble creating vertices and quads for Euler grid.");
-}
-
-ErrorCode NCHelperEuler::init_EulSpcscd_vals(const FileOptions& opts, EntityHandle file_set)
-{
Interface*& mbImpl = _readNC->mbImpl;
std::vector<std::string>& dimNames = _readNC->dimNames;
std::vector<int>& dimVals = _readNC->dimVals;
@@ -472,7 +473,7 @@ ErrorCode NCHelperEuler::init_EulSpcscd_vals(const FileOptions& opts, EntityHand
return MB_SUCCESS;
}
-ErrorCode NCHelperEuler::create_EulSpcscd_verts_quads(ScdInterface* scdi, EntityHandle file_set, Range& quads)
+ErrorCode NCHelperEuler::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
{
return _readNC->create_scd_verts_quads(scdi, file_set, quads);
}
diff --git a/src/io/NCHelperEuler.hpp b/src/io/NCHelperEuler.hpp
index affea4a..b1c4c21 100644
--- a/src/io/NCHelperEuler.hpp
+++ b/src/io/NCHelperEuler.hpp
@@ -19,7 +19,7 @@ class NCHelperEuler : public NCHelper
public:
NCHelperEuler(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
- static bool can_read_file(ReadNC* readNC);
+ static bool can_read_file(ReadNC* readNC, int fileId);
private:
virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
@@ -29,9 +29,6 @@ private:
virtual std::string get_mesh_type_name() { return "CAM_EUL"; }
virtual bool is_scd_mesh() { return true; }
-
- ErrorCode init_EulSpcscd_vals(const FileOptions& opts, EntityHandle file_set);
- ErrorCode create_EulSpcscd_verts_quads(ScdInterface* scdi, EntityHandle file_set, Range& quads);
};
} // namespace moab
diff --git a/src/io/NCHelperFV.cpp b/src/io/NCHelperFV.cpp
index ff42fbe..b086102 100644
--- a/src/io/NCHelperFV.cpp
+++ b/src/io/NCHelperFV.cpp
@@ -10,39 +10,40 @@
namespace moab {
-bool NCHelperFV::can_read_file(ReadNC* readNC)
+bool NCHelperFV::can_read_file(ReadNC* readNC, int fileId)
{
std::vector<std::string>& dimNames = readNC->dimNames;
- // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it's the FV grid
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it should be the FV grid
if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
dimNames.end(), std::string("lat")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slon"))
- != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end()))
+ != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end())) {
+ // Make sure it is CAM grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
+ if (attIt == readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
+ return false;
+ }
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ if (success != 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
+ return false;
+ }
+ if (att_data.find("CAM") == std::string::npos)
+ return false;
+
return true;
+ }
return false;
}
ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
{
- ErrorCode rval = init_FVCDscd_vals(opts, file_set);
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble initializing FV grid.");
-
- return rval;
-}
-
-ErrorCode NCHelperFV::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
-{
- ErrorCode rval = create_FVCDscd_verts_quads(scdi, file_set, quads);
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble creating vertices and quads for FV grid.");
-
- return rval;
-}
-
-ErrorCode NCHelperFV::init_FVCDscd_vals(const FileOptions& opts, EntityHandle file_set)
-{
Interface*& mbImpl = _readNC->mbImpl;
std::vector<std::string>& dimNames = _readNC->dimNames;
std::vector<int>& dimVals = _readNC->dimVals;
@@ -121,7 +122,7 @@ ErrorCode NCHelperFV::init_FVCDscd_vals(const FileOptions& opts, EntityHandle fi
if (globallyPeriodic[0])
assert("Number of vertices and edges should be same" && gDims[3] == gCDims[3]);
else
- assert("Number of vertices should equal to number of edges plus one" && gDims[3] == gCDims[3]+1);
+ assert("Number of vertices should equal to number of edges plus one" && gDims[3] == gCDims[3] + 1);
#ifdef USE_MPI
// if serial, use a locally-periodic representation only if local mesh is periodic, otherwise don't
@@ -468,7 +469,7 @@ ErrorCode NCHelperFV::init_FVCDscd_vals(const FileOptions& opts, EntityHandle fi
return MB_SUCCESS;
}
-ErrorCode NCHelperFV::create_FVCDscd_verts_quads(ScdInterface* scdi, EntityHandle file_set, Range& quads)
+ErrorCode NCHelperFV::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
{
return _readNC->create_scd_verts_quads(scdi, file_set, quads);
}
diff --git a/src/io/NCHelperFV.hpp b/src/io/NCHelperFV.hpp
index b9f5d31..afcd9b8 100644
--- a/src/io/NCHelperFV.hpp
+++ b/src/io/NCHelperFV.hpp
@@ -18,17 +18,13 @@ class NCHelperFV : public NCHelper
{
public:
NCHelperFV(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
- static bool can_read_file(ReadNC* readNC);
+ static bool can_read_file(ReadNC* readNC, int fileId);
private:
virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads);
virtual std::string get_mesh_type_name() { return "CAM_FV"; }
virtual bool is_scd_mesh() { return true; }
-
-private:
- ErrorCode init_FVCDscd_vals(const FileOptions& opts, EntityHandle file_set);
- ErrorCode create_FVCDscd_verts_quads(ScdInterface* scdi, EntityHandle file_set, Range& quads);
};
} // namespace moab
diff --git a/src/io/NCHelperHOMME.cpp b/src/io/NCHelperHOMME.cpp
index 604e000..6490b24 100644
--- a/src/io/NCHelperHOMME.cpp
+++ b/src/io/NCHelperHOMME.cpp
@@ -13,48 +13,51 @@
namespace moab {
-bool NCHelperHOMME::can_read_file(ReadNC* readNC, const FileOptions& opts, int& spectralOrder)
+NCHelperHOMME::NCHelperHOMME(ReadNC* readNC, int fileId, const FileOptions& opts) : NCHelper(readNC, fileId), _spectralOrder(-1)
{
- // If global attribute "np" exists then it's the HOMME grid
+ // Calculate spectral order
std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("np");
- if (attIt != readNC->globalAtts.end())
- {
- int success = NCFUNC(get_att_int)(readNC->fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &spectralOrder);
- if (success != 0) {
+ if (attIt != readNC->globalAtts.end()) {
+ int success = NCFUNC(get_att_int)(readNC->fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &_spectralOrder);
+ if (success != 0)
readNC->readMeshIface->report_error("%s", "Failed to read np global attribute int data.");
- return false;
- }
-
- spectralOrder--; // Spectral order is one less than np
+ else
+ _spectralOrder--; // Spectral order is one less than np
if (MB_SUCCESS == opts.match_option("PARTITION_METHOD", "NODAL_PARTITION"))
readNC->partMethod = -1;
-
- return true;
}
-
- return false;
}
-ErrorCode NCHelperHOMME::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+bool NCHelperHOMME::can_read_file(ReadNC* readNC, int fileId)
{
- ErrorCode rval = init_HOMMEucd_vals();
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble initializing HOMME grid.");
-
- return rval;
-}
+ // If global attribute "np" exists then it should be the HOMME grid
+ if (readNC->globalAtts.find("np") != readNC->globalAtts.end()) {
+ // Make sure it is CAM grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
+ if (attIt == readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
+ return false;
+ }
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ if (success != 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
+ return false;
+ }
+ if (att_data.find("CAM") == std::string::npos)
+ return false;
-ErrorCode NCHelperHOMME::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
-{
- ErrorCode rval = create_HOMMEucd_verts_quads(opts, file_set, quads);
- if (MB_SUCCESS != rval)
- _readNC->readMeshIface->report_error("%s", "Trouble creating vertices and quads for HOMME grid.");
+ return true;
+ }
- return rval;
+ return false;
}
-ErrorCode NCHelperHOMME::init_HOMMEucd_vals()
+ErrorCode NCHelperHOMME::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
{
std::vector<std::string>& dimNames = _readNC->dimNames;
std::vector<int>& dimVals = _readNC->dimVals;
@@ -80,7 +83,9 @@ ErrorCode NCHelperHOMME::init_HOMMEucd_vals()
idx = vit - dimNames.begin();
else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
idx = vit - dimNames.begin();
- else ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ }
tDim = idx;
tMax = dimVals[idx] - 1;
tMin = 0;
@@ -186,10 +191,15 @@ ErrorCode NCHelperHOMME::init_HOMMEucd_vals()
// with no corresponding variables
_readNC->init_dims_with_no_cvars_info();
+ // This check is for HOMME and other ucd mesh. When ReadNC class instance
+ // gets out of scope in a script (and deleted), the localGid will be lost
+ rval = _readNC->check_ucd_localGid(file_set);
+ ERRORR(rval, "Trouble checking local Gid for ucd mesh.");
+
return MB_SUCCESS;
}
-ErrorCode NCHelperHOMME::create_HOMMEucd_verts_quads(const FileOptions& opts, EntityHandle file_set, Range& quads)
+ErrorCode NCHelperHOMME::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
{
Interface*& mbImpl = _readNC->mbImpl;
std::string& fileName = _readNC->fileName;
@@ -260,15 +270,10 @@ ErrorCode NCHelperHOMME::create_HOMMEucd_verts_quads(const FileOptions& opts, En
std::vector<std::string>::iterator vit;
int idx;
if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncells")) != conn_names.end())
- {
idx = vit - conn_names.begin();
- }
else if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncenters")) != conn_names.end())
- {
idx = vit - conn_names.begin();
- }
- else
- {
+ else {
ERRORR(MB_FAILURE, "Failed to get number of quads.");
}
int num_quads = conn_vals[idx];
@@ -277,7 +282,7 @@ ErrorCode NCHelperHOMME::create_HOMMEucd_verts_quads(const FileOptions& opts, En
int cornerVarId;
success = NCFUNC(inq_varid)(connectId, "element_corners", &cornerVarId);
ERRORS(success, "Failed to get variable id.");
- NCDF_SIZE tmp_dims[2] = { 0, 0 }, tmp_counts[2] = { 4, static_cast<size_t>(num_quads) };
+ NCDF_SIZE tmp_dims[2] = {0, 0}, tmp_counts[2] = {4, static_cast<size_t>(num_quads)};
std::vector<int> tmp_conn(4 * num_quads), tmp_conn2(4 * num_quads);
success = NCFUNCAG(_vara_int)(connectId, cornerVarId, tmp_dims, tmp_counts, &tmp_conn2[0] NCREQ);
ERRORS(success, "Failed to get temporary connectivity.");
@@ -285,10 +290,10 @@ ErrorCode NCHelperHOMME::create_HOMMEucd_verts_quads(const FileOptions& opts, En
ERRORS(success, "Failed on close.");
// permute the connectivity
for (int i = 0; i < num_quads; i++) {
- tmp_conn[4*i] = tmp_conn2[i];
- tmp_conn[4*i + 1] = tmp_conn2[i + 1 * num_quads];
- tmp_conn[4*i + 2] = tmp_conn2[i + 2 * num_quads];
- tmp_conn[4*i + 3] = tmp_conn2[i + 3 * num_quads];
+ tmp_conn[4 * i] = tmp_conn2[i];
+ tmp_conn[4 * i + 1] = tmp_conn2[i + 1 * num_quads];
+ tmp_conn[4 * i + 2] = tmp_conn2[i + 2 * num_quads];
+ tmp_conn[4 * i + 3] = tmp_conn2[i + 3 * num_quads];
}
// need to know whether we'll be creating gather mesh later, to make sure we allocate enough space
@@ -309,7 +314,8 @@ ErrorCode NCHelperHOMME::create_HOMMEucd_verts_quads(const FileOptions& opts, En
start_idx = 4 * rank * num_coarse_quads * spectral_unit;
// iextra = # coarse quads extra after equal split over procs
int iextra = num_quads % (procs * spectral_unit);
- if (rank < iextra) num_coarse_quads++;
+ if (rank < iextra)
+ num_coarse_quads++;
start_idx += 4 * spectral_unit * std::min(rank, iextra);
// num_fine_quads is the number of quads in the connectivity list in HommeMapping file assigned to this proc
num_fine_quads = spectral_unit * num_coarse_quads;
@@ -359,7 +365,7 @@ ErrorCode NCHelperHOMME::create_HOMMEucd_verts_quads(const FileOptions& opts, En
double *xptr = arrays[0], *yptr = arrays[1], *zptr = arrays[2];
int i;
for (i = 0, rit = localGid.begin(); i < (int)num_local_verts; i++, ++rit) {
- assert(*rit < ilVals.size()+1);
+ assert(*rit < ilVals.size() + 1);
xptr[i] = ilVals[(*rit) - 1];
yptr[i] = jlVals[(*rit) - 1];
zptr[i] = klVals[lDims[2]];
@@ -383,8 +389,7 @@ ErrorCode NCHelperHOMME::create_HOMMEucd_verts_quads(const FileOptions& opts, En
int *gid_data = (int*) data;
std::copy(localGid.begin(), localGid.end(), gid_data);
// duplicate global id data, which will be used to resolve sharing
- if (mpFileIdTag)
- {
+ if (mpFileIdTag) {
rval = mbImpl->tag_iterate(*mpFileIdTag, vert_range.begin(), vert_range.end(), count, data);
ERRORR(rval, "Failed to get tag iterator on file id tag.");
assert(count == (int) num_local_verts);
diff --git a/src/io/NCHelperHOMME.hpp b/src/io/NCHelperHOMME.hpp
index f0dc1cd..feb3cad 100644
--- a/src/io/NCHelperHOMME.hpp
+++ b/src/io/NCHelperHOMME.hpp
@@ -17,8 +17,8 @@ namespace moab {
class NCHelperHOMME : public NCHelper
{
public:
- NCHelperHOMME(ReadNC* readNC, int fileId, int spectralOrder) : NCHelper(readNC, fileId), _spectralOrder(spectralOrder) {}
- static bool can_read_file(ReadNC* readNC, const FileOptions& opts, int& spectralOrder);
+ NCHelperHOMME(ReadNC* readNC, int fileId, const FileOptions& opts);
+ static bool can_read_file(ReadNC* readNC, int fileId);
private:
virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
@@ -26,9 +26,6 @@ private:
virtual std::string get_mesh_type_name() { return "CAM_SE"; }
virtual bool is_scd_mesh() { return false; }
- ErrorCode init_HOMMEucd_vals();
- ErrorCode create_HOMMEucd_verts_quads(const FileOptions& opts, EntityHandle file_set, Range& quads);
-
private:
int _spectralOrder; // read from variable 'np'
};
diff --git a/src/io/ReadNC.cpp b/src/io/ReadNC.cpp
index 89e0494..8c77431 100644
--- a/src/io/ReadNC.cpp
+++ b/src/io/ReadNC.cpp
@@ -203,8 +203,7 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
}
rval = helper->init_mesh_vals(opts, tmp_set);
- if (MB_SUCCESS != rval)
- return rval;
+ ERRORR(rval, "Trouble initializing mesh values.");
// Create mesh vertex/quads sequences
Range quads;
@@ -214,33 +213,10 @@ ErrorCode ReadNC::load_file(const char *file_name, const EntityHandle* file_set,
}
else if (!noMesh) {
rval = helper->create_verts_quads(scdi, opts, tmp_set, quads);
- if (MB_SUCCESS != rval)
- return rval;
+ ERRORR(rval, "Trouble creating vertices and quads.");
}
bool scdMesh = helper->is_scd_mesh();
- if (noMesh && !scdMesh)
- {
- // we need to populate localGid range with the gids of vertices from the tmp_set
- // localGid is important in reading the variable data into the nodes
- // also, for our purposes, localGid is truly the GLOBAL_ID tag data, not other
- // file_id tags that could get passed around in other scenarios for parallel reading
- // for nodal_partition, this local gid is easier, should be initialized with only
- // the owned nodes
-
- // we need to get all vertices from tmp_set (it is the input set in no_mesh scenario)
- Range local_verts;
- rval = mbImpl->get_entities_by_dimension(tmp_set, 0, local_verts);
- if (MB_FAILURE == rval)
- return rval;
- std::vector<int> gids(local_verts.size());
- // !IMPORTANT : this has to be the GLOBAL_ID tag
- rval=mbImpl->tag_get_data(mGlobalIdTag, local_verts, &gids[0]);
- if (MB_FAILURE == rval)
- return rval;
- // this will do a smart copy
- std::copy(gids.begin(), gids.end(), range_inserter(localGid));
- }
// Read variables onto grid
if (!noVars) {
@@ -465,6 +441,41 @@ ErrorCode ReadNC::parse_options(const FileOptions &opts, std::vector<std::string
return MB_SUCCESS;
}
+// In a script, the ReadNC class instance can get out of scope (and deleted). In that
+// case, the localGid (initialized properly when the mesh was created) will be lost,
+// so it has to be properly refilled with the Global Ids of the local vertices
+ErrorCode ReadNC::check_ucd_localGid(EntityHandle tmp_set)
+{
+ if (noMesh && localGid.empty()) {
+ // we need to populate localGid range with the gids of vertices from the tmp_set
+ // localGid is important in reading the variable data into the nodes
+ // also, for our purposes, localGid is truly the GLOBAL_ID tag data, not other
+ // file_id tags that could get passed around in other scenarios for parallel reading
+ // for nodal_partition, this local gid is easier, should be initialized with only
+ // the owned nodes
+
+ // we need to get all vertices from tmp_set (it is the input set in no_mesh scenario)
+ Range local_verts;
+ ErrorCode rval = mbImpl->get_entities_by_dimension(tmp_set, 0, local_verts);
+ if (MB_FAILURE == rval)
+ return rval;
+
+ if (!local_verts.empty()) {
+ std::vector<int> gids(local_verts.size());
+
+ // !IMPORTANT : this has to be the GLOBAL_ID tag
+ rval = mbImpl->tag_get_data(mGlobalIdTag, local_verts, &gids[0]);
+ if (MB_FAILURE == rval)
+ return rval;
+
+ // this will do a smart copy
+ std::copy(gids.begin(), gids.end(), range_inserter(localGid));
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
ErrorCode ReadNC::check_verts_quads(EntityHandle file_set) {
// check parameters on this read against what was on the mesh from last read
// get the number of vertices
diff --git a/src/io/ReadNC.hpp b/src/io/ReadNC.hpp
index 61d4224..87ed049 100644
--- a/src/io/ReadNC.hpp
+++ b/src/io/ReadNC.hpp
@@ -161,6 +161,9 @@ private:
//! create vertices for scd mesh
ErrorCode create_scd_verts_quads(ScdInterface *scdi, EntityHandle file_set, Range &quads);
+ //! make sure that localGid is properly initialized for ucd mesh
+ ErrorCode check_ucd_localGid(EntityHandle file_set);
+
//! check number of vertices and elements against what's already in file_set
ErrorCode check_verts_quads(EntityHandle file_set);
https://bitbucket.org/fathomteam/moab/commits/4a602615bf50/
Changeset: 4a602615bf50
Branch: master
User: tautges
Date: 2013-06-05 18:49:22
Summary: Merged in danwu/moab (pull request #5)
Update ReadNC and NCHelper based on last code review
Affected #: 11 files
diff --git a/src/io/Makefile.am b/src/io/Makefile.am
index 9c76524..b7eae6a 100644
--- a/src/io/Makefile.am
+++ b/src/io/Makefile.am
@@ -16,8 +16,12 @@ if NETCDF_FILE
MOAB_NETCDF_SRCS = ReadNCDF.cpp ReadNCDF.hpp \
WriteNCDF.cpp WriteNCDF.hpp \
WriteSLAC.cpp WriteSLAC.hpp \
- ReadNC.cpp ReadNC.hpp \
- ReadGCRM.cpp ReadGCRM.hpp
+ ReadNC.cpp ReadNC.hpp \
+ NCHelper.cpp NCHelper.hpp \
+ NCHelperEuler.cpp NCHelperEuler.hpp \
+ NCHelperFV.cpp NCHelperFV.hpp \
+ NCHelperHOMME.cpp NCHelperHOMME.hpp \
+ ReadGCRM.cpp ReadGCRM.hpp
else
MOAB_NETCDF_SRCS =
endif
@@ -25,7 +29,11 @@ endif
if PNETCDF_FILE
if !NETCDF_FILE
MOAB_NETCDF_SRCS += ReadNC.cpp ReadNC.hpp \
- ReadGCRM.cpp ReadGCRM.hpp
+ NCHelper.cpp NCHelper.hpp \
+ NCHelperEuler.cpp NCHelperEuler.hpp \
+ NCHelperFV.cpp NCHelperFV.hpp \
+ NCHelperHOMME.cpp NCHelperHOMME.hpp \
+ ReadGCRM.cpp ReadGCRM.hpp
endif
endif
diff --git a/src/io/NCHelper.cpp b/src/io/NCHelper.cpp
new file mode 100644
index 0000000..e0f021e
--- /dev/null
+++ b/src/io/NCHelper.cpp
@@ -0,0 +1,21 @@
+#include "NCHelper.hpp"
+#include "NCHelperEuler.hpp"
+#include "NCHelperFV.hpp"
+#include "NCHelperHOMME.hpp"
+#include "moab/ReadUtilIface.hpp"
+
+namespace moab {
+
+NCHelper* NCHelper::get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts)
+{
+ if (NCHelperEuler::can_read_file(readNC, fileId))
+ return new (std::nothrow) NCHelperEuler(readNC, fileId);
+ else if (NCHelperFV::can_read_file(readNC, fileId))
+ return new (std::nothrow) NCHelperFV(readNC, fileId);
+ else if (NCHelperHOMME::can_read_file(readNC, fileId))
+ return new (std::nothrow) NCHelperHOMME(readNC, fileId, opts);
+ else // Unknown NetCDF grid (will fill this in later for POP, CICE and CLM)
+ return NULL;
+}
+
+} // namespace moab
diff --git a/src/io/NCHelper.hpp b/src/io/NCHelper.hpp
new file mode 100644
index 0000000..852fc61
--- /dev/null
+++ b/src/io/NCHelper.hpp
@@ -0,0 +1,36 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelper.hpp
+//
+// Purpose : Climate NC file helper
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPER_HPP
+#define NCHELPER_HPP
+
+#include "ReadNC.hpp"
+
+namespace moab {
+
+//! Helper class to isolate reading of several different nc file formats
+class NCHelper
+{
+public:
+ NCHelper(ReadNC* readNC, int fileId) : _readNC(readNC), _fileId(fileId) {}
+
+ static NCHelper* get_nc_helper(ReadNC* readNC, int fileId, const FileOptions& opts);
+
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set) = 0;
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads) = 0;
+ virtual std::string get_mesh_type_name() = 0;
+ virtual bool is_scd_mesh() = 0;
+
+protected:
+ ReadNC* _readNC;
+ int _fileId;
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCHelperEuler.cpp b/src/io/NCHelperEuler.cpp
new file mode 100644
index 0000000..9346d22
--- /dev/null
+++ b/src/io/NCHelperEuler.cpp
@@ -0,0 +1,481 @@
+#include "NCHelperEuler.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+
+#include <cmath>
+#include <sstream>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) {_readNC->readMeshIface->report_error("%s", str); return rval;}
+
+#define ERRORS(err, str) \
+ if (err) {_readNC->readMeshIface->report_error("%s", str); return MB_FAILURE;}
+
+namespace moab {
+
+bool NCHelperEuler::can_read_file(ReadNC* readNC, int fileId)
+{
+ std::vector<std::string>& dimNames = readNC->dimNames;
+
+ // If dimension names "lon" AND "lat' exist then it could be either the Eulerian Spectral grid or the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) != dimNames.end())) {
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it should be the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("slon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("slat")) != dimNames.end()))
+ return false;
+
+ // Make sure it is CAM grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
+ if (attIt == readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
+ return false;
+ }
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ if (success != 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
+ return false;
+ }
+ if (att_data.find("CAM") == std::string::npos)
+ return false;
+
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ Interface*& mbImpl = _readNC->mbImpl;
+ std::vector<std::string>& dimNames = _readNC->dimNames;
+ std::vector<int>& dimVals = _readNC->dimVals;
+ std::string& iName = _readNC->iName;
+ std::string& jName = _readNC->jName;
+ std::string& tName = _readNC->tName;
+ std::string& iCName = _readNC->iCName;
+ std::string& jCName = _readNC->jCName;
+ std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
+ int& tMin = _readNC->tMin;
+ int& tMax = _readNC->tMax;
+ int (&gDims)[6] = _readNC->gDims;
+ int (&lDims)[6] = _readNC->lDims;
+ int (&lCDims)[6] = _readNC->lCDims;
+ int (&gCDims)[6] = _readNC->gCDims;
+ std::vector<double>& ilVals = _readNC->ilVals;
+ std::vector<double>& jlVals = _readNC->jlVals;
+ std::vector<double>& tVals = _readNC->tVals;
+ std::vector<double>& ilCVals = _readNC->ilCVals;
+ std::vector<double>& jlCVals = _readNC->jlCVals;
+ int& tDim = _readNC->tDim;
+ int& iCDim = _readNC->iCDim;
+ int& jCDim = _readNC->jCDim;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+ bool& isParallel = _readNC->isParallel;
+ int& partMethod = _readNC->partMethod;
+ int (&locallyPeriodic)[2] = _readNC->locallyPeriodic;
+ int (&globallyPeriodic)[2] = _readNC->globallyPeriodic;
+ ScdParData& parData = _readNC->parData;
+#ifdef USE_MPI
+ ParallelComm*& myPcomm = _readNC->myPcomm;
+#endif
+
+ // look for names of center i/j dimensions
+ std::vector<std::string>::iterator vit;
+ unsigned int idx;
+ iCName = std::string("lon");
+ iName = std::string("slon");
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), iCName.c_str())) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find center i variable.");
+ }
+ iCDim = idx;
+
+ // decide on i periodicity using math for now
+ std::vector<double> tilVals(dimVals[idx]);
+ ErrorCode rval = _readNC->read_coordinate(iCName.c_str(), 0, dimVals[idx] - 1, tilVals);
+ ERRORR(rval, "Trouble reading lon variable.");
+ if (std::fabs(2 * (*(tilVals.rbegin())) - *(tilVals.rbegin() + 1) - 360) < 0.001)
+ globallyPeriodic[0] = 1;
+
+ // now we can set gCDims and gDims for i
+ gCDims[0] = 0;
+ gDims[0] = 0;
+ gCDims[3] = dimVals[idx] - 1; // these are stored directly in file
+ gDims[3] = gCDims[3] + (globallyPeriodic[0] ? 0 : 1); // only if not periodic is vertex param max > elem param max
+
+ // now j
+ jCName = std::string("lat");
+ jName = std::string("slat");
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), jCName.c_str())) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find center j variable.");
+ }
+ jCDim = idx;
+
+ // for Eul models, will always be non-periodic in j
+ gCDims[1] = 0;
+ gDims[1] = 0;
+ gCDims[4] = dimVals[idx] - 1;
+ gDims[4] = gCDims[4] + 1;
+
+ // try a truly 2d mesh
+ gDims[2] = -1;
+ gDims[5] = -1;
+
+ // look for time dimensions
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ }
+ tDim = idx;
+ tMax = dimVals[idx] - 1;
+ tMin = 0;
+ tName = dimNames[idx];
+
+ // parse options to get subset
+ if (isParallel) {
+#ifdef USE_MPI
+ for (int i = 0; i < 6; i++)
+ parData.gDims[i] = gDims[i];
+ for (int i = 0; i < 3; i++)
+ parData.gPeriodic[i] = globallyPeriodic[i];
+ parData.partMethod = partMethod;
+ int pdims[3];
+
+ rval = ScdInterface::compute_partition(myPcomm->proc_config().proc_size(),
+ myPcomm->proc_config().proc_rank(),
+ parData, lDims, locallyPeriodic, pdims);
+ if (MB_SUCCESS != rval)
+ return rval;
+ for (int i = 0; i < 3; i++)
+ parData.pDims[i] = pdims[i];
+
+ dbgOut.tprintf(1, "Partition: %dx%d (out of %dx%d)\n",
+ lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
+ gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1);
+ if (myPcomm->proc_config().proc_rank() == 0)
+ dbgOut.tprintf(1, "Contiguous chunks of size %d bytes.\n", 8 * (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1));
+#endif
+ }
+ else {
+ for (int i = 0; i < 6; i++)
+ lDims[i] = gDims[i];
+ locallyPeriodic[0] = globallyPeriodic[0];
+ }
+
+ opts.get_int_option("IMIN", lDims[0]);
+ opts.get_int_option("IMAX", lDims[3]);
+ opts.get_int_option("JMIN", lDims[1]);
+ opts.get_int_option("JMAX", lDims[4]);
+
+ // now get actual coordinate values for vertices and cell centers; first resize
+ if (locallyPeriodic[0]) {
+ // if locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem coords
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0] + 1);
+ lCDims[3] = lDims[3];
+ }
+ else {
+ if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3]) {
+ // globally periodic and I'm the last proc, get fewer vertex coords than vertices in i
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] = lDims[3] - 1;
+ }
+ else {
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] = lDims[3] - 1;
+ }
+ }
+
+ lCDims[0] = lDims[0];
+ lCDims[4] = lDims[4] - 1;
+ lCDims[1] = lDims[1];
+
+ if (-1 != lDims[1]) {
+ jlVals.resize(lDims[4] - lDims[1] + 1);
+ jlCVals.resize(lCDims[4] - lCDims[1] + 1);
+ }
+
+ if (-1 != tMin)
+ tVals.resize(tMax - tMin + 1);
+
+ // now read coord values
+ std::map<std::string, ReadNC::VarData>::iterator vmit;
+ if (!ilCVals.empty()) {
+ if ((vmit = varInfo.find(iCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(iCName.c_str(), lDims[0], lDims[0] + ilCVals.size() - 1, ilCVals);
+ ERRORR(rval, "Trouble reading lon variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lon coordinate.");
+ }
+ }
+
+ if (!jlCVals.empty()) {
+ if ((vmit = varInfo.find(jCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(jCName.c_str(), lDims[1], lDims[1] + jlCVals.size() - 1, jlCVals);
+ ERRORR(rval, "Trouble reading lat variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lat coordinate.");
+ }
+ }
+
+ if (lDims[0] != -1) {
+ if ((vmit = varInfo.find(iCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ double dif = (ilCVals[1] - ilCVals[0]) / 2;
+ std::size_t i;
+ for (i = 0; i != ilCVals.size(); i++)
+ ilVals[i] = ilCVals[i] - dif;
+ // the last one is needed only if not periodic
+ if (!locallyPeriodic[0])
+ ilVals[i] = ilCVals[i - 1] + dif;
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ }
+ }
+
+ if (lDims[1] != -1) {
+ if ((vmit = varInfo.find(jCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ if (!isParallel || ((gDims[4] - gDims[1]) == (lDims[4] - lDims[1]))) {
+ std::string gwName("gw");
+ std::vector<double> gwVals(lDims[4] - lDims[1] - 1);
+ rval = _readNC->read_coordinate(gwName.c_str(), lDims[1], lDims[4] - 2, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ // copy the correct piece
+ jlVals[0] = -(M_PI / 2) * 180 / M_PI;
+ unsigned int i = 0;
+ double gwSum = -1;
+ for (i = 1; i != gwVals.size() + 1; i++) {
+ gwSum = gwSum + gwVals[i - 1];
+ jlVals[i] = std::asin(gwSum) * 180 / M_PI;
+ }
+ jlVals[i] = 90.0; // using value of i after loop exits.
+ }
+ else {
+ std::string gwName("gw");
+ double gwSum = 0;
+
+ // If this is the first row
+ if (lDims[1] == gDims[1]) {
+ std::vector<double> gwVals(lDims[4]);
+ rval = _readNC->read_coordinate(gwName.c_str(), 0, lDims[4] - 1, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ // copy the correct piece
+ jlVals[0] = -(M_PI / 2) * 180 / M_PI;
+ gwSum = -1;
+ for (std::size_t i = 1; i != jlVals.size(); i++) {
+ gwSum = gwSum + gwVals[i - 1];
+ jlVals[i] = std::asin(gwSum) * 180 / M_PI;
+ }
+ }
+ // or if it's the last row
+ else if (lDims[4] == gDims[4]) {
+ std::vector<double> gwVals(lDims[4] - 1);
+ rval = _readNC->read_coordinate(gwName.c_str(), 0, lDims[4] - 2, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ // copy the correct piece
+ gwSum = -1;
+ for (int j = 0; j != lDims[1] - 1; j++) {
+ gwSum = gwSum + gwVals[j];
+ }
+ std::size_t i = 0;
+ for (; i != jlVals.size() - 1; i++) {
+ gwSum = gwSum + gwVals[lDims[1] - 1 + i];
+ jlVals[i] = std::asin(gwSum) * 180 / M_PI;
+ }
+ jlVals[i] = 90.0; // using value of i after loop exits.
+ }
+ // it's in the middle
+ else {
+ int start = lDims[1] - 1;
+ int end = lDims[4] - 1;
+ std::vector<double> gwVals(end);
+ rval = _readNC->read_coordinate(gwName.c_str(), 0, end - 1, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ gwSum = -1;
+ for (int j = 0; j != start - 1; j++) {
+ gwSum = gwSum + gwVals[j];
+ }
+ std::size_t i = 0;
+ for (; i != jlVals.size(); i++) {
+ gwSum = gwSum + gwVals[start - 1 + i];
+ jlVals[i] = std::asin(gwSum) * 180 / M_PI;
+ }
+ }
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ }
+ }
+
+ if (tMin != -1) {
+ if ((vmit = varInfo.find(tName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+ }
+
+ dbgOut.tprintf(1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4]);
+ dbgOut.tprintf(1, "%d elements, %d vertices\n", (lDims[3] - lDims[0]) * (lDims[4] - lDims[1]), (lDims[3] - lDims[0] + 1)
+ * (lDims[4] - lDims[1] + 1));
+
+ // determine the entity location type of a variable
+ std::map<std::string, ReadNC::VarData>::iterator mit;
+ for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) {
+ ReadNC::VarData& vd = (*mit).second;
+ if ((std::find(vd.varDims.begin(), vd.varDims.end(), iCDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), jCDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCQUAD;
+ }
+
+ // <coordinate_dim_name>
+ std::vector<std::string> ijdimNames(4);
+ ijdimNames[0] = "__slon";
+ ijdimNames[1] = "__slat";
+ ijdimNames[2] = "__lon";
+ ijdimNames[3] = "__lat";
+
+ std::string tag_name;
+ int val_len = 0;
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ tag_name = ijdimNames[i];
+ void * val = NULL;
+ if (tag_name == "__slon") {
+ val = &ilVals[0];
+ val_len = ilVals.size();
+ }
+ else if (tag_name == "__slat") {
+ val = &jlVals[0];
+ val_len = jlVals.size();
+ }
+ else if (tag_name == "__lon") {
+ val = &ilCVals[0];
+ val_len = ilCVals.size();
+ }
+ else if (tag_name == "__lat") {
+ val = &jlCVals[0];
+ val_len = jlCVals.size();
+ }
+ Tag tagh = 0;
+ DataType data_type;
+
+ // assume all has same data type as lon
+ switch (varInfo["lon"].varDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ case NC_DOUBLE:
+ data_type = MB_TYPE_DOUBLE;
+ break;
+ case NC_FLOAT:
+ data_type = MB_TYPE_DOUBLE;
+ break;
+ case NC_INT:
+ data_type = MB_TYPE_INTEGER;
+ break;
+ case NC_SHORT:
+ default:
+ std::cerr << "Unrecognized data type for tag " << tag_name << std::endl;
+ ERRORR(MB_FAILURE, "Unrecognized data type");
+ break;
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, data_type, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating <coordinate_dim_name> tag.");
+ rval = mbImpl->tag_set_by_ptr(tagh, &file_set, 1, &val, &val_len);
+ ERRORR(rval, "Trouble setting data for <coordinate_dim_name> tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_LOC_MINMAX
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
+ tag_name = ss_tag_name.str();
+ Tag tagh = 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] == "__slon") {
+ val[0] = lDims[0];
+ val[1] = lDims[3];
+ }
+ else if (ijdimNames[i] == "__slat") {
+ val[0] = lDims[1];
+ val[1] = lDims[4];
+ }
+ else if (ijdimNames[i] == "__lon") {
+ val[0] = lCDims[0];
+ val[1] = lCDims[3];
+ }
+ else if (ijdimNames[i] == "__lat") {
+ val[0] = lCDims[1];
+ val[1] = lCDims[4];
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_LOC_MINMAX tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_LOC_MINMAX tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_GLOBAL_MINMAX
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
+ tag_name = ss_tag_name.str();
+ Tag tagh = 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] == "__slon") {
+ val[0] = gDims[0];
+ val[1] = gDims[3];
+ }
+ else if (ijdimNames[i] == "__slat") {
+ val[0] = gDims[1];
+ val[1] = gDims[4];
+ }
+ else if (ijdimNames[i] == "__lon") {
+ val[0] = gCDims[0];
+ val[1] = gCDims[3];
+ }
+ else if (ijdimNames[i] == "__lat") {
+ val[0] = gCDims[1];
+ val[1] = gCDims[4];
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_GLOBAL_MINMAX tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_GLOBAL_MINMAX tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // hack: create dummy tags, if needed, for variables like nbnd
+ // with no corresponding variables
+ _readNC->init_dims_with_no_cvars_info();
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelperEuler::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+{
+ return _readNC->create_scd_verts_quads(scdi, file_set, quads);
+}
+
+} // namespace moab
diff --git a/src/io/NCHelperEuler.hpp b/src/io/NCHelperEuler.hpp
new file mode 100644
index 0000000..b1c4c21
--- /dev/null
+++ b/src/io/NCHelperEuler.hpp
@@ -0,0 +1,36 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelperEuler.hpp
+//
+// Purpose : Climate NC file helper for Eulerian Spectral grid
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPEREULER_HPP
+#define NCHELPEREULER_HPP
+
+#include "NCHelper.hpp"
+
+namespace moab {
+
+//! Child helper class for Eulerian Spectral grid (CAM_EUL)
+class NCHelperEuler : public NCHelper
+{
+public:
+ NCHelperEuler(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+ static bool can_read_file(ReadNC* readNC, int fileId);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads);
+
+ virtual std::string get_mesh_type_name() { return "CAM_EUL"; }
+
+ virtual bool is_scd_mesh() { return true; }
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCHelperFV.cpp b/src/io/NCHelperFV.cpp
new file mode 100644
index 0000000..b086102
--- /dev/null
+++ b/src/io/NCHelperFV.cpp
@@ -0,0 +1,477 @@
+#include "NCHelperFV.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+
+#include <cmath>
+#include <sstream>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) {_readNC->readMeshIface->report_error("%s", str); return rval;}
+
+namespace moab {
+
+bool NCHelperFV::can_read_file(ReadNC* readNC, int fileId)
+{
+ std::vector<std::string>& dimNames = readNC->dimNames;
+
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it should be the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) != dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slon"))
+ != dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(), std::string("slat")) != dimNames.end())) {
+ // Make sure it is CAM grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
+ if (attIt == readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
+ return false;
+ }
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ if (success != 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
+ return false;
+ }
+ if (att_data.find("CAM") == std::string::npos)
+ return false;
+
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ Interface*& mbImpl = _readNC->mbImpl;
+ std::vector<std::string>& dimNames = _readNC->dimNames;
+ std::vector<int>& dimVals = _readNC->dimVals;
+ std::string& iName = _readNC->iName;
+ std::string& jName = _readNC->jName;
+ std::string& tName = _readNC->tName;
+ std::string& iCName = _readNC->iCName;
+ std::string& jCName = _readNC->jCName;
+ std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
+ int& tMin = _readNC->tMin;
+ int& tMax = _readNC->tMax;
+ int (&gDims)[6] = _readNC->gDims;
+ int (&lDims)[6] = _readNC->lDims;
+ int (&gCDims)[6] = _readNC->gCDims;
+ int (&lCDims)[6] = _readNC->lCDims;
+ std::vector<double>& ilVals = _readNC->ilVals;
+ std::vector<double>& jlVals = _readNC->jlVals;
+ std::vector<double>& tVals = _readNC->tVals;
+ std::vector<double>& ilCVals = _readNC->ilCVals;
+ std::vector<double>& jlCVals = _readNC->jlCVals;
+ int& iDim = _readNC->iDim;
+ int& jDim = _readNC->jDim;
+ int& tDim = _readNC->tDim;
+ int& iCDim = _readNC->iCDim;
+ int& jCDim = _readNC->jCDim;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+ bool& isParallel = _readNC->isParallel;
+ int& partMethod = _readNC->partMethod;
+ int (&locallyPeriodic)[2] = _readNC->locallyPeriodic;
+ int (&globallyPeriodic)[2] = _readNC->globallyPeriodic;
+ ScdParData& parData = _readNC->parData;
+#ifdef USE_MPI
+ ParallelComm*& myPcomm = _readNC->myPcomm;
+#endif
+
+ std::vector<std::string>::iterator vit;
+ unsigned int idx;
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "slon")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find slon variable.");
+ }
+
+ iDim = idx;
+ gDims[3] = dimVals[idx] - 1;
+ gDims[0] = 0;
+ iName = dimNames[idx];
+
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "slat")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find slat variable.");
+ }
+ jDim = idx;
+ gDims[4] = dimVals[idx] - 1 + 2; // add 2 for the pole points
+ gDims[1] = 0;
+ jName = dimNames[idx];
+
+ // look for names of center i/j dimensions
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lon")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lon variable.");
+ }
+ iCDim = idx;
+ gCDims[3] = dimVals[idx] - 1;
+ gCDims[0] = 0;
+ iCName = dimNames[idx];
+
+ // check and set globallyPeriodic[0]
+ std::vector<double> til_vals(2);
+ ErrorCode rval = _readNC->read_coordinate(iCName.c_str(), dimVals[idx] - 2, dimVals[idx] - 1, til_vals);
+ ERRORR(rval, "Trouble reading slon variable.");
+ if (std::fabs(2 * til_vals[1] - til_vals[0] - 360) < 0.001)
+ globallyPeriodic[0] = 1;
+ if (globallyPeriodic[0])
+ assert("Number of vertices and edges should be same" && gDims[3] == gCDims[3]);
+ else
+ assert("Number of vertices should equal to number of edges plus one" && gDims[3] == gCDims[3] + 1);
+
+#ifdef USE_MPI
+ // if serial, use a locally-periodic representation only if local mesh is periodic, otherwise don't
+ if ((isParallel && myPcomm->proc_config().proc_size() == 1) && globallyPeriodic[0])
+ locallyPeriodic[0] = 1;
+#else
+ if (globallyPeriodic[0])
+ locallyPeriodic[0] = 1;
+#endif
+
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lat")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lat variable.");
+ }
+ jCDim = idx;
+ gCDims[4] = dimVals[idx] - 1;
+ gCDims[1] = 0;
+ jCName = dimNames[idx];
+
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ }
+ tDim = idx;
+ tMax = dimVals[idx] - 1;
+ tMin = 0;
+ tName = dimNames[idx];
+
+ // parse options to get subset
+ if (isParallel) {
+#ifdef USE_MPI
+ for (int i = 0; i < 6; i++)
+ parData.gDims[i] = gDims[i];
+ for (int i = 0; i < 3; i++)
+ parData.gPeriodic[i] = globallyPeriodic[i];
+ parData.partMethod = partMethod;
+ int pdims[3];
+
+ rval = ScdInterface::compute_partition(myPcomm->proc_config().proc_size(),
+ myPcomm->proc_config().proc_rank(),
+ parData, lDims, locallyPeriodic, pdims);
+ if (MB_SUCCESS != rval)
+ return rval;
+ for (int i = 0; i < 3; i++)
+ parData.pDims[i] = pdims[i];
+
+ dbgOut.tprintf(1, "Partition: %dx%d (out of %dx%d)\n",
+ lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
+ gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1);
+ if (myPcomm->proc_config().proc_rank() == 0)
+ dbgOut.tprintf(1, "Contiguous chunks of size %d bytes.\n", 8 * (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1));
+#endif
+ }
+ else {
+ for (int i = 0; i < 6; i++)
+ lDims[i] = gDims[i];
+ locallyPeriodic[0] = globallyPeriodic[0];
+ }
+ opts.get_int_option("IMIN", lDims[0]);
+ opts.get_int_option("IMAX", lDims[3]);
+ opts.get_int_option("JMIN", lDims[1]);
+ opts.get_int_option("JMAX", lDims[4]);
+
+ // now get actual coordinate values for vertices and cell centers; first resize
+ if (locallyPeriodic[0]) {
+ // if locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem coords
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0] + 1);
+ lCDims[3] = lDims[3];
+ }
+ else {
+ if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3]) {
+ // globally periodic and I'm the last proc, get fewer vertex coords than vertices in i
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] = lDims[3] - 1;
+ }
+ else {
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] = lDims[3] - 1;
+ }
+ }
+
+ lCDims[0] = lDims[0];
+ lCDims[4] = lDims[4] - 1;
+ lCDims[1] = lDims[1];
+
+ if (-1 != lDims[1]) {
+ jlVals.resize(lDims[4] - lDims[1] + 1);
+ jlCVals.resize(lCDims[4] - lCDims[1] + 1);
+ }
+
+ if (-1 != tMin)
+ tVals.resize(tMax - tMin + 1);
+
+ // ... then read actual values
+ std::map<std::string, ReadNC::VarData>::iterator vmit;
+ if (lCDims[0] != -1) {
+ if ((vmit = varInfo.find(iCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(iCName.c_str(), lCDims[0], lCDims[3], ilCVals);
+ ERRORR(rval, "Trouble reading lon variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lon coordinate.");
+ }
+ }
+
+ if (lCDims[1] != -1) {
+ if ((vmit = varInfo.find(jCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(jCName.c_str(), lCDims[1], lCDims[4], jlCVals);
+ ERRORR(rval, "Trouble reading lat variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lat coordinate.");
+ }
+ }
+
+ if (lDims[0] != -1) {
+ if ((vmit = varInfo.find(iName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ // last column
+ if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3]) {
+ til_vals.resize(ilVals.size() - 1, 0.0);
+ rval = _readNC->read_coordinate(iName.c_str(), lDims[0], lDims[3] - 1, til_vals);
+ double dif = til_vals[1] - til_vals[0];
+ std::size_t i;
+ for (i = 0; i != til_vals.size(); i++)
+ ilVals[i] = til_vals[i];
+ ilVals[i] = ilVals[i - 1] + dif;
+ }
+ else {
+ rval = _readNC->read_coordinate(iName.c_str(), lDims[0], lDims[3], ilVals);
+ ERRORR(rval, "Trouble reading x variable.");
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ }
+ }
+
+ if (lDims[1] != -1) {
+ if ((vmit = varInfo.find(jName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ if (!isParallel || ((gDims[4] - gDims[1]) == (lDims[4] - lDims[1]))) {
+ std::vector<double> dummyVar(lDims[4] - lDims[1] - 1);
+ rval = _readNC->read_coordinate(jName.c_str(), lDims[1], lDims[4] - 2, dummyVar);
+ ERRORR(rval, "Trouble reading y variable.");
+ // copy the correct piece
+ jlVals[0] = -90.0;
+ unsigned int i = 0;
+ for (i = 1; i != dummyVar.size() + 1; i++)
+ jlVals[i] = dummyVar[i - 1];
+ jlVals[i] = 90.0; // using value of i after loop exits.
+ }
+ else {
+ // If this is the first row
+ // need to read one less then available and read it into a dummy var
+ if (lDims[1] == gDims[1]) {
+ std::vector<double> dummyVar(lDims[4] - lDims[1]);
+ rval = _readNC->read_coordinate(jName.c_str(), lDims[1], lDims[4] - 1, dummyVar);
+ ERRORR(rval, "Trouble reading y variable.");
+ // copy the correct piece
+ jlVals[0] = -90.0;
+ for (int i = 1; i < lDims[4] + 1; i++)
+ jlVals[i] = dummyVar[i - 1];
+ }
+ // or if it's the last row
+ else if (lDims[4] == gDims[4]) {
+ std::vector<double> dummyVar(lDims[4] - lDims[1]);
+ rval = _readNC->read_coordinate(jName.c_str(), lDims[1] - 1, lDims[4] - 2, dummyVar);
+ ERRORR(rval, "Trouble reading y variable.");
+ // copy the correct piece
+ std::size_t i = 0;
+ for (i = 0; i != dummyVar.size(); i++)
+ jlVals[i] = dummyVar[i];
+ jlVals[i] = 90.0; // using value of i after loop exits.
+ }
+ // it's in the middle
+ else {
+ rval = _readNC->read_coordinate(jCName.c_str(), lDims[1] - 1, lDims[4] - 1, jlVals);
+ ERRORR(rval, "Trouble reading y variable.");
+ }
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ }
+ }
+
+ if (tMin != -1) {
+ if ((vmit = varInfo.find(tName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+ }
+
+ dbgOut.tprintf(1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4]);
+ dbgOut.tprintf(1, "%d elements, %d vertices\n", (lDims[3] - lDims[0]) * (lDims[4] - lDims[1]), (lDims[3] - lDims[0] + 1)
+ * (lDims[4] - lDims[1] + 1));
+
+ // determine the entity location type of a variable
+ std::map<std::string, ReadNC::VarData>::iterator mit;
+ for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) {
+ ReadNC::VarData& vd = (*mit).second;
+ if ((std::find(vd.varDims.begin(), vd.varDims.end(), iCDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), jCDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCQUAD;
+ else if ((std::find(vd.varDims.begin(), vd.varDims.end(), jDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), iCDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCNSEDGE;
+ else if ((std::find(vd.varDims.begin(), vd.varDims.end(), jCDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), iDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCEWEDGE;
+ }
+
+ // <coordinate_dim_name>
+ std::vector<std::string> ijdimNames(4);
+ ijdimNames[0] = "__slon";
+ ijdimNames[1] = "__slat";
+ ijdimNames[2] = "__lon";
+ ijdimNames[3] = "__lat";
+
+ std::string tag_name;
+ int val_len = 0;
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ tag_name = ijdimNames[i];
+ void * val = NULL;
+ if (tag_name == "__slon") {
+ val = &ilVals[0];
+ val_len = ilVals.size();
+ }
+ else if (tag_name == "__slat") {
+ val = &jlVals[0];
+ val_len = jlVals.size();
+ }
+ else if (tag_name == "__lon") {
+ val = &ilCVals[0];
+ val_len = ilCVals.size();
+ }
+ else if (tag_name == "__lat") {
+ val = &jlCVals[0];
+ val_len = jlCVals.size();
+ }
+ Tag tagh = 0;
+ DataType data_type;
+
+ // assume all has same data type as lon
+ switch (varInfo["lon"].varDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ case NC_DOUBLE:
+ data_type = MB_TYPE_DOUBLE;
+ break;
+ case NC_FLOAT:
+ data_type = MB_TYPE_DOUBLE;
+ break;
+ case NC_INT:
+ data_type = MB_TYPE_INTEGER;
+ break;
+ case NC_SHORT:
+ default:
+ std::cerr << "Unrecognized data type for tag " << tag_name << std::endl;
+ ERRORR(MB_FAILURE, "Unrecognized data type");
+ break;
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, data_type, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating <coordinate_dim_name> tag.");
+ rval = mbImpl->tag_set_by_ptr(tagh, &file_set, 1, &val, &val_len);
+ ERRORR(rval, "Trouble setting data for <coordinate_dim_name> tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_LOC_MINMAX
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
+ tag_name = ss_tag_name.str();
+ Tag tagh = 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] == "__slon") {
+ val[0] = lDims[0];
+ val[1] = lDims[3];
+ }
+ else if (ijdimNames[i] == "__slat") {
+ val[0] = lDims[1];
+ val[1] = lDims[4];
+ }
+ else if (ijdimNames[i] == "__lon") {
+ val[0] = lCDims[0];
+ val[1] = lCDims[3];
+ }
+ else if (ijdimNames[i] == "__lat") {
+ val[0] = lCDims[1];
+ val[1] = lCDims[4];
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_LOC_MINMAX tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_LOC_MINMAX tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_GLOBAL_MINMAX
+ for (unsigned int i = 0; i != ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
+ tag_name = ss_tag_name.str();
+ Tag tagh = 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] == "__slon") {
+ val[0] = gDims[0];
+ val[1] = gDims[3];
+ }
+ else if (ijdimNames[i] == "__slat") {
+ val[0] = gDims[1];
+ val[1] = gDims[4];
+ }
+ else if (ijdimNames[i] == "__lon") {
+ val[0] = gCDims[0];
+ val[1] = gCDims[3];
+ }
+ else if (ijdimNames[i] == "__lat") {
+ val[0] = gCDims[1];
+ val[1] = gCDims[4];
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_GLOBAL_MINMAX tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_GLOBAL_MINMAX tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // hack: create dummy tags, if needed, for dimensions like nbnd
+ // with no corresponding variables
+ _readNC->init_dims_with_no_cvars_info();
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelperFV::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+{
+ return _readNC->create_scd_verts_quads(scdi, file_set, quads);
+}
+
+} // namespace moab
diff --git a/src/io/NCHelperFV.hpp b/src/io/NCHelperFV.hpp
new file mode 100644
index 0000000..afcd9b8
--- /dev/null
+++ b/src/io/NCHelperFV.hpp
@@ -0,0 +1,32 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelperFV.hpp
+//
+// Purpose : Climate NC file helper for Finite Volume grid
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPERFV_HPP
+#define NCHELPERFV_HPP
+
+#include "NCHelper.hpp"
+
+namespace moab {
+
+//! Child helper class for Finite Volume grid (CAM_FV)
+class NCHelperFV : public NCHelper
+{
+public:
+ NCHelperFV(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+ static bool can_read_file(ReadNC* readNC, int fileId);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads);
+ virtual std::string get_mesh_type_name() { return "CAM_FV"; }
+ virtual bool is_scd_mesh() { return true; }
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCHelperHOMME.cpp b/src/io/NCHelperHOMME.cpp
new file mode 100644
index 0000000..6490b24
--- /dev/null
+++ b/src/io/NCHelperHOMME.cpp
@@ -0,0 +1,506 @@
+#include "NCHelperHOMME.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+#include "moab/SpectralMeshTool.hpp"
+
+#include <cmath>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) {_readNC->readMeshIface->report_error("%s", str); return rval;}
+
+#define ERRORS(err, str) \
+ if (err) {_readNC->readMeshIface->report_error("%s", str); return MB_FAILURE;}
+
+namespace moab {
+
+NCHelperHOMME::NCHelperHOMME(ReadNC* readNC, int fileId, const FileOptions& opts) : NCHelper(readNC, fileId), _spectralOrder(-1)
+{
+ // Calculate spectral order
+ std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("np");
+ if (attIt != readNC->globalAtts.end()) {
+ int success = NCFUNC(get_att_int)(readNC->fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &_spectralOrder);
+ if (success != 0)
+ readNC->readMeshIface->report_error("%s", "Failed to read np global attribute int data.");
+ else
+ _spectralOrder--; // Spectral order is one less than np
+
+ if (MB_SUCCESS == opts.match_option("PARTITION_METHOD", "NODAL_PARTITION"))
+ readNC->partMethod = -1;
+ }
+}
+
+bool NCHelperHOMME::can_read_file(ReadNC* readNC, int fileId)
+{
+ // If global attribute "np" exists then it should be the HOMME grid
+ if (readNC->globalAtts.find("np") != readNC->globalAtts.end()) {
+ // Make sure it is CAM grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
+ if (attIt == readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
+ return false;
+ }
+ unsigned int sz = attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] = '\000';
+ int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
+ if (success != 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
+ return false;
+ }
+ if (att_data.find("CAM") == std::string::npos)
+ return false;
+
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHelperHOMME::init_mesh_vals(const FileOptions& opts, EntityHandle file_set)
+{
+ std::vector<std::string>& dimNames = _readNC->dimNames;
+ std::vector<int>& dimVals = _readNC->dimVals;
+ std::string& kName = _readNC->kName;
+ std::string& tName = _readNC->tName;
+ std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
+ int& tMin = _readNC->tMin;
+ int& tMax = _readNC->tMax;
+ int (&gDims)[6] = _readNC->gDims;
+ int (&lDims)[6] = _readNC->lDims;
+ int& iDim = _readNC->iDim;
+ int& kDim = _readNC->kDim;
+ int& tDim = _readNC->tDim;
+ std::vector<double>& ilVals = _readNC->ilVals;
+ std::vector<double>& jlVals = _readNC->jlVals;
+ std::vector<double>& klVals = _readNC->klVals;
+ std::vector<double>& tVals = _readNC->tVals;
+
+ ErrorCode rval;
+ unsigned int idx;
+ std::vector<std::string>::iterator vit;
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ }
+ tDim = idx;
+ tMax = dimVals[idx] - 1;
+ tMin = 0;
+ tName = dimNames[idx];
+
+ // get number of vertices (labeled as number of columns) and levels
+ gDims[0] = gDims[3] = -1;
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "ncol")) != dimNames.end()) {
+ idx = vit - dimNames.begin();
+ gDims[3] = dimVals[idx] - 1;
+ gDims[0] = 0;
+ iDim = idx;
+ }
+ if (-1 == gDims[0])
+ return MB_FAILURE;
+
+ // set j coordinate to the number of quads
+ gDims[1] = gDims[0];
+ gDims[4] = gDims[3] - 2;
+
+ gDims[2] = gDims[5] = -1;
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end()) {
+ idx = vit - dimNames.begin();
+ gDims[5] = dimVals[idx] - 1, gDims[2] = 0, kName = std::string("lev");
+ kDim = idx;
+ }
+ if (-1 == gDims[2])
+ return MB_FAILURE;
+
+ // read coordinate data
+ std::map<std::string, ReadNC::VarData>::iterator vmit;
+ if (gDims[0] != -1) {
+ if ((vmit = varInfo.find("lon")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate("lon", gDims[0], gDims[3], ilVals);
+ ERRORR(rval, "Trouble reading x variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ }
+ }
+
+ // store lat values in jlVals parameterized by j
+ if (gDims[1] != -1) {
+ if ((vmit = varInfo.find("lat")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate("lat", gDims[0], gDims[3], jlVals);
+ ERRORR(rval, "Trouble reading y variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ }
+ }
+
+ if (gDims[2] != -1) {
+ if ((vmit = varInfo.find("lev")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate("lev", gDims[2], gDims[5], klVals);
+ ERRORR(rval, "Trouble reading z variable.");
+
+ // decide whether down is positive
+ char posval[10];
+ int success = NCFUNC(get_att_text)(_readNC->fileId, (*vmit).second.varId, "positive", posval);
+ if (0 == success && !strcmp(posval, "down")) {
+ for (std::vector<double>::iterator dvit = klVals.begin(); dvit != klVals.end(); ++dvit)
+ (*dvit) *= -1.0;
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find z coordinate.");
+ }
+ }
+
+ if (tMin != -1) {
+ if ((vmit = varInfo.find(tName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+ }
+
+ if ((vmit = varInfo.find(tName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+
+ // determine the entity location type of a variable
+ std::map<std::string, ReadNC::VarData>::iterator mit;
+ for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) {
+ ReadNC::VarData& vd = (*mit).second;
+ if ((std::find(vd.varDims.begin(), vd.varDims.end(), iDim) != vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), kDim) != vd.varDims.end()))
+ vd.entLoc = ReadNC::ENTLOCNODE;
+ }
+
+ std::copy(gDims, gDims + 6, lDims);
+
+ // don't read coordinates of columns until we actually create the mesh
+
+ // hack: create dummy tags, if needed, for variables like ncol and nbnd
+ // with no corresponding variables
+ _readNC->init_dims_with_no_cvars_info();
+
+ // This check is for HOMME and other ucd mesh. When ReadNC class instance
+ // gets out of scope in a script (and deleted), the localGid will be lost
+ rval = _readNC->check_ucd_localGid(file_set);
+ ERRORR(rval, "Trouble checking local Gid for ucd mesh.");
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelperHOMME::create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads)
+{
+ Interface*& mbImpl = _readNC->mbImpl;
+ std::string& fileName = _readNC->fileName;
+ int& connectId = _readNC->connectId;
+ int (&gDims)[6] = _readNC->gDims;
+ int (&lDims)[6] = _readNC->lDims;
+ std::vector<double>& ilVals = _readNC->ilVals;
+ std::vector<double>& jlVals = _readNC->jlVals;
+ std::vector<double>& klVals = _readNC->klVals;
+ Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
+ const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+ bool& isParallel = _readNC->isParallel;
+ Range& localGid = _readNC->localGid;
+#ifdef USE_MPI
+ ParallelComm*& myPcomm = _readNC->myPcomm;
+#endif
+ bool& spectralMesh = _readNC->spectralMesh;
+
+ // need to get/read connectivity data before creating elements
+ std::string conn_fname;
+
+ // try to open the connectivity file through CONN option, if used
+ ErrorCode rval = opts.get_str_option("CONN", conn_fname);
+ if (MB_SUCCESS != rval) {
+ // default convention for reading HOMME is a file HommeMapping.nc in same dir as data file
+ conn_fname = std::string(fileName);
+ size_t idx = conn_fname.find_last_of("/");
+ if (idx != std::string::npos)
+ conn_fname = conn_fname.substr(0, idx).append("/HommeMapping.nc");
+ else
+ conn_fname = "HommeMapping.nc";
+ }
+
+ int success;
+
+ int rank, procs;
+#ifdef PNETCDF_FILE
+ if (isParallel) {
+ success = NCFUNC(open)(myPcomm->proc_config().proc_comm(), conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId);
+ rank = myPcomm->proc_config().proc_rank();
+ procs = myPcomm->proc_config().proc_size();
+ }
+ else {
+ success = NCFUNC(open)(MPI_COMM_SELF, conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId);
+ rank = 0;
+ procs = 1;
+ }
+#else
+ success = NCFUNC(open)(conn_fname.c_str(), 0, &connectId);
+ rank = 0;
+ procs = 1;
+#endif
+ ERRORS(success, "Failed on open.");
+
+ std::vector<std::string> conn_names;
+ std::vector<int> conn_vals;
+ rval = _readNC->get_dimensions(connectId, conn_names, conn_vals);
+ ERRORR(rval, "Failed to get dimensions for connectivity.");
+
+ if (conn_vals[0] != gDims[3] - gDims[0] + 1 - 2) {
+ dbgOut.tprintf(1, "Warning: number of quads from %s and vertices from %s are inconsistent; nverts = %d, nquads = %d.\n",
+ conn_fname.c_str(), fileName.c_str(), gDims[3] - gDims[0] + 1, conn_vals[0]);
+ }
+
+ // read connectivity into temporary variable
+ int num_fine_quads, num_coarse_quads, start_idx;
+ std::vector<std::string>::iterator vit;
+ int idx;
+ if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncells")) != conn_names.end())
+ idx = vit - conn_names.begin();
+ else if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncenters")) != conn_names.end())
+ idx = vit - conn_names.begin();
+ else {
+ ERRORR(MB_FAILURE, "Failed to get number of quads.");
+ }
+ int num_quads = conn_vals[idx];
+
+ // get the connectivity into tmp_conn2 and permute into tmp_conn
+ int cornerVarId;
+ success = NCFUNC(inq_varid)(connectId, "element_corners", &cornerVarId);
+ ERRORS(success, "Failed to get variable id.");
+ NCDF_SIZE tmp_dims[2] = {0, 0}, tmp_counts[2] = {4, static_cast<size_t>(num_quads)};
+ std::vector<int> tmp_conn(4 * num_quads), tmp_conn2(4 * num_quads);
+ success = NCFUNCAG(_vara_int)(connectId, cornerVarId, tmp_dims, tmp_counts, &tmp_conn2[0] NCREQ);
+ ERRORS(success, "Failed to get temporary connectivity.");
+ success = NCFUNC(close)(connectId);
+ ERRORS(success, "Failed on close.");
+ // permute the connectivity
+ for (int i = 0; i < num_quads; i++) {
+ tmp_conn[4 * i] = tmp_conn2[i];
+ tmp_conn[4 * i + 1] = tmp_conn2[i + 1 * num_quads];
+ tmp_conn[4 * i + 2] = tmp_conn2[i + 2 * num_quads];
+ tmp_conn[4 * i + 3] = tmp_conn2[i + 3 * num_quads];
+ }
+
+ // need to know whether we'll be creating gather mesh later, to make sure we allocate enough space
+ // in one shot
+ bool create_gathers = true;
+#ifdef USE_MPI
+ if (isParallel)
+ if (myPcomm->proc_config().proc_rank() != 0)
+ create_gathers = false;
+#endif
+
+ // compute the number of local quads, accounting for coarse or fine representation
+ // spectral_unit is the # fine quads per coarse quad, or spectralOrder^2
+ int spectral_unit = (spectralMesh ? _spectralOrder * _spectralOrder : 1);
+ // num_coarse_quads is the number of quads instantiated in MOAB; if !spectralMesh, num_coarse_quads = num_fine_quads
+ num_coarse_quads = int(std::floor(1.0 * num_quads / (spectral_unit * procs)));
+ // start_idx is the starting index in the HommeMapping connectivity list for this proc, before converting to coarse quad representation
+ start_idx = 4 * rank * num_coarse_quads * spectral_unit;
+ // iextra = # coarse quads extra after equal split over procs
+ int iextra = num_quads % (procs * spectral_unit);
+ if (rank < iextra)
+ num_coarse_quads++;
+ start_idx += 4 * spectral_unit * std::min(rank, iextra);
+ // num_fine_quads is the number of quads in the connectivity list in HommeMapping file assigned to this proc
+ num_fine_quads = spectral_unit * num_coarse_quads;
+
+ // now create num_coarse_quads
+ EntityHandle *conn_arr;
+ EntityHandle start_vertex;
+ Range tmp_range;
+
+ // read connectivity into that space
+ EntityHandle *sv_ptr = NULL, start_quad;
+ SpectralMeshTool smt(mbImpl, _spectralOrder);
+ if (!spectralMesh) {
+ rval = _readNC->readMeshIface->get_element_connect(num_coarse_quads, 4,
+ MBQUAD, 0, start_quad, conn_arr,
+ // might have to create gather mesh later
+ (create_gathers ? num_coarse_quads + num_quads : num_coarse_quads));
+ ERRORR(rval, "Failed to create quads.");
+ tmp_range.insert(start_quad, start_quad + num_coarse_quads - 1);
+ std::copy(&tmp_conn[start_idx], &tmp_conn[start_idx + 4 * num_fine_quads], conn_arr);
+ std::copy(conn_arr, conn_arr + 4 * num_fine_quads, range_inserter(localGid));
+ }
+ else {
+ rval = smt.create_spectral_elems(&tmp_conn[0], num_fine_quads, 2, tmp_range, start_idx, &localGid);
+ ERRORR(rval, "Failed to create spectral elements.");
+ int count, v_per_e;
+ rval = mbImpl->connect_iterate(tmp_range.begin(), tmp_range.end(), conn_arr, v_per_e, count);
+ ERRORR(rval, "Failed to get connectivity of spectral elements.");
+ rval = mbImpl->tag_iterate(smt.spectral_vertices_tag(true), tmp_range.begin(), tmp_range.end(),
+ count, (void*&)sv_ptr);
+ ERRORR(rval, "Failed to get fine connectivity of spectral elements.");
+ }
+
+ // on this proc, I get columns ldims[1]..ldims[4], inclusive; need to find which vertices those correpond to
+ unsigned int num_local_verts = localGid.size();
+ unsigned int num_total_verts = gDims[3] - gDims[0] + 1;
+
+ // create vertices
+ std::vector<double*> arrays;
+ rval = _readNC->readMeshIface->get_node_coords(3, num_local_verts, 0, start_vertex, arrays,
+ // might have to create gather mesh later
+ (create_gathers ? num_local_verts+num_total_verts : num_local_verts));
+ ERRORR(rval, "Couldn't create vertices in ucd mesh.");
+
+ // set vertex coordinates
+ Range::iterator rit;
+ double *xptr = arrays[0], *yptr = arrays[1], *zptr = arrays[2];
+ int i;
+ for (i = 0, rit = localGid.begin(); i < (int)num_local_verts; i++, ++rit) {
+ assert(*rit < ilVals.size() + 1);
+ xptr[i] = ilVals[(*rit) - 1];
+ yptr[i] = jlVals[(*rit) - 1];
+ zptr[i] = klVals[lDims[2]];
+ }
+
+ const double pideg = acos(-1.0) / 180.0;
+ for (i = 0; i < (int)num_local_verts; i++) {
+ double cosphi = cos(pideg * yptr[i]);
+ double zmult = sin(pideg * yptr[i]), xmult = cosphi * cos(xptr[i] * pideg), ymult = cosphi * sin(xptr[i] * pideg);
+ double rad = 8.0e3 + klVals[lDims[2]];
+ xptr[i] = rad * xmult, yptr[i] = rad * ymult, zptr[i] = rad * zmult;
+ }
+
+ // get ptr to gid memory for vertices
+ Range vert_range(start_vertex, start_vertex + num_local_verts - 1);
+ void *data;
+ int count;
+ rval = mbImpl->tag_iterate(mGlobalIdTag, vert_range.begin(), vert_range.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator.");
+ assert(count == (int) num_local_verts);
+ int *gid_data = (int*) data;
+ std::copy(localGid.begin(), localGid.end(), gid_data);
+ // duplicate global id data, which will be used to resolve sharing
+ if (mpFileIdTag) {
+ rval = mbImpl->tag_iterate(*mpFileIdTag, vert_range.begin(), vert_range.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator on file id tag.");
+ assert(count == (int) num_local_verts);
+ gid_data = (int*) data;
+ std::copy(localGid.begin(), localGid.end(), gid_data);
+ }
+
+ // create map from file ids to vertex handles, used later to set connectivity
+ std::map<EntityHandle, EntityHandle> vert_handles;
+ for (rit = localGid.begin(), i = 0; rit != localGid.end(); ++rit, i++) {
+ vert_handles[*rit] = start_vertex + i;
+ }
+
+ // compute proper handles in connectivity using offset
+ for (int q = 0; q < 4 * num_coarse_quads; q++) {
+ conn_arr[q] = vert_handles[conn_arr[q]];
+ assert(conn_arr[q]);
+ }
+ if (spectralMesh) {
+ int verts_per_quad = (_spectralOrder + 1) * (_spectralOrder + 1);
+ for (int q = 0; q < verts_per_quad * num_coarse_quads; q++) {
+ sv_ptr[q] = vert_handles[sv_ptr[q]];
+ assert(sv_ptr[q]);
+ }
+ }
+
+ // add new vertices and elements to the set
+ quads.merge(tmp_range);
+ tmp_range.insert(start_vertex, start_vertex + num_local_verts - 1);
+ rval = mbImpl->add_entities(file_set, tmp_range);
+ ERRORR(rval, "Couldn't add new vertices and quads/hexes to file set.");
+
+ // mark the set with the spectral order
+ Tag sporder;
+ rval = mbImpl->tag_get_handle("SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, sporder, MB_TAG_CREAT | MB_TAG_SPARSE);
+ ERRORR(rval, "Couldn't create spectral order tag.");
+ rval = mbImpl->tag_set_data(sporder, &file_set, 1, &_spectralOrder);
+ ERRORR(rval, "Couldn't set value for spectral order tag.");
+
+#ifdef USE_MPI
+ if (isParallel && myPcomm->proc_config().proc_rank() == 0) {
+#endif
+ EntityHandle gather_set;
+ rval = mbImpl->create_meshset(MESHSET_SET, gather_set);
+ ERRORR(rval, "Trouble creating gather set.");
+
+ // create vertices
+ arrays.clear();
+ // don't need to specify allocation number here, because we know enough verts were created before
+ rval = _readNC->readMeshIface->get_node_coords(3, num_total_verts, 0, start_vertex, arrays);
+ ERRORR(rval, "Couldn't create vertices in ucd mesh for gather set.");
+
+ xptr = arrays[0], yptr = arrays[1], zptr = arrays[2];
+ for (i = 0; i < (int)num_total_verts; i++) {
+ double cosphi = cos(pideg * jlVals[i]);
+ double zmult = sin(pideg * jlVals[i]);
+ double xmult = cosphi * cos(ilVals[i] * pideg);
+ double ymult = cosphi * sin(ilVals[i] * pideg);
+ double rad = 8.0e3 + klVals[lDims[2]];
+ xptr[i] = rad * xmult;
+ yptr[i] = rad * ymult;
+ zptr[i] = rad * zmult;
+ }
+
+ // get ptr to gid memory for vertices
+ Range gather_verts(start_vertex, start_vertex + num_total_verts - 1);
+ rval = mbImpl->tag_iterate(mGlobalIdTag, gather_verts.begin(), gather_verts.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator.");
+ assert(count == (int) num_total_verts);
+ gid_data = (int*) data;
+ for (int j = 1; j <= (int) num_total_verts; j++)
+ gid_data[j - 1] = j;
+ // set the file id tag too, it should be bigger something not interfering with global id
+ if (mpFileIdTag) {
+ rval = mbImpl->tag_iterate(*mpFileIdTag, gather_verts.begin(), gather_verts.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator in file id tag.");
+ assert(count == (int) num_total_verts);
+ gid_data = (int*) data;
+ for (int j = 1; j <= (int) num_total_verts; j++)
+ gid_data[j - 1] = num_total_verts + j; // bigger than global id tag
+ }
+
+ rval = mbImpl->add_entities(gather_set, gather_verts);
+ ERRORR(rval, "Couldn't add vertices to gather set.");
+
+ // create quads
+ Range gather_quads;
+ // don't need to specify allocation number here, because we know enough quads were created before
+ rval = _readNC->readMeshIface->get_element_connect(num_quads, 4,
+ MBQUAD, 0, start_quad, conn_arr);
+ ERRORR(rval, "Failed to create quads.");
+ gather_quads.insert(start_quad, start_quad + num_quads - 1);
+ std::copy(&tmp_conn[0], &tmp_conn[4 * num_quads], conn_arr);
+ for (i = 0; i != 4 * num_quads; i++)
+ conn_arr[i] += start_vertex - 1; // connectivity array is shifted by where the gather verts start
+ rval = mbImpl->add_entities(gather_set, gather_quads);
+ ERRORR(rval, "Couldn't add quads to gather set.");
+
+ Tag gathersettag;
+ rval = mbImpl->tag_get_handle("GATHER_SET", 1, MB_TYPE_INTEGER, gathersettag,
+ MB_TAG_CREAT | MB_TAG_SPARSE);
+ ERRORR(rval, "Couldn't create gather set tag.");
+ int gatherval = 1;
+ rval = mbImpl->tag_set_data(gathersettag, &gather_set, 1, &gatherval);
+ ERRORR(rval, "Couldn't set value for gather set tag.");
+
+#ifdef USE_MPI
+ }
+#endif
+
+ return MB_SUCCESS;
+}
+
+} // namespace moab
diff --git a/src/io/NCHelperHOMME.hpp b/src/io/NCHelperHOMME.hpp
new file mode 100644
index 0000000..feb3cad
--- /dev/null
+++ b/src/io/NCHelperHOMME.hpp
@@ -0,0 +1,35 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelperHOMME.hpp
+//
+// Purpose : Climate NC file helper for HOMME grid
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPERHOMME_HPP
+#define NCHELPERHOMME_HPP
+
+#include "NCHelper.hpp"
+
+namespace moab {
+
+//! Child helper class for HOMME grid (CAM_SE)
+class NCHelperHOMME : public NCHelper
+{
+public:
+ NCHelperHOMME(ReadNC* readNC, int fileId, const FileOptions& opts);
+ static bool can_read_file(ReadNC* readNC, int fileId);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle file_set);
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& quads);
+ virtual std::string get_mesh_type_name() { return "CAM_SE"; }
+ virtual bool is_scd_mesh() { return false; }
+
+private:
+ int _spectralOrder; // read from variable 'np'
+};
+
+} // namespace moab
+
+#endif
This diff is so big that we needed to truncate the remainder.
Repository URL: https://bitbucket.org/fathomteam/moab/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the moab-dev
mailing list