[MOAB-dev] commit/MOAB: 2 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Aug 12 14:53:35 CDT 2013
2 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/d011d558c33a/
Changeset: d011d558c33a
Branch: None
User: tautges
Date: 2013-08-12 20:52:45
Summary: Changes to ScdInterface, mostly to support new partitioning method (sqijk) but also to allow
for periodic meshes in all 3 directions.
NOTE: all applications using structured mesh functions in parallel should check parameters to
ScdInterface::compute_partition and ScdInterface::get_neighbor; the parameters for periodicity and
across_bdy have changed dimension from 2 to 3.
Also, fixed a problem in scdseq_test where returned value didn't reflect #errors, so tests
were passing when they should have been failing.
Affected #: 10 files
diff --git a/src/ScdInterface.cpp b/src/ScdInterface.cpp
index 8084640..7c91ab7 100644
--- a/src/ScdInterface.cpp
+++ b/src/ScdInterface.cpp
@@ -248,7 +248,7 @@ Tag ScdInterface::box_periodic_tag(bool create_if_missing)
{
if (boxPeriodicTag || !create_if_missing) return boxPeriodicTag;
- ErrorCode rval = mbImpl->tag_get_handle("BOX_PERIODIC", 2, MB_TYPE_INTEGER,
+ ErrorCode rval = mbImpl->tag_get_handle("BOX_PERIODIC", 3, MB_TYPE_INTEGER,
boxPeriodicTag, MB_TAG_SPARSE|MB_TAG_CREAT);
if (MB_SUCCESS != rval) return 0;
return boxPeriodicTag;
@@ -693,7 +693,7 @@ ErrorCode ScdInterface::get_neighbor_alljkbal(int np, int pfrom,
}
pto = -1;
- across_bdy[0] = across_bdy[1] = 0;
+ across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
int ldims[6], pijk[3], lperiodic[2];
ErrorCode rval = compute_partition_alljkbal(np, pfrom, gdims, gperiodic,
@@ -801,7 +801,7 @@ ErrorCode ScdInterface::get_neighbor_sqij(int np, int pfrom,
}
pto = -1;
- across_bdy[0] = across_bdy[1] = 0;
+ across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
int lperiodic[3], pijk[3], ldims[6];
ErrorCode rval = compute_partition_sqij(np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk);
if (MB_SUCCESS != rval) return rval;
@@ -929,7 +929,7 @@ ErrorCode ScdInterface::get_neighbor_sqjk(int np, int pfrom,
}
pto = -1;
- across_bdy[0] = across_bdy[1] = 0;
+ across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
int pijk[3], lperiodic[3], ldims[6];
ErrorCode rval = compute_partition_sqjk(np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk);
if (MB_SUCCESS != rval) return rval;
@@ -952,7 +952,7 @@ ErrorCode ScdInterface::get_neighbor_sqjk(int np, int pfrom,
pto = pfrom;
int dj = (gdims[4] - gdims[1]) / pijk[1], jextra = (gdims[4] - gdims[1]) % dj,
dk = (gdims[5] == gdims[2] ? 0 : (gdims[5] - gdims[2]) / pijk[2]), kextra = (gdims[5] - gdims[2]) - dk*pijk[2];
-
+ assert((dj*pijk[1] + jextra == (gdims[4]-gdims[1])) && (dk*pijk[2] + kextra == (gdims[5]-gdims[2])));
if (0 != dijk[1]) {
pto = (nj + dijk[1] + pijk[1]) % pijk[1]; // get pto's ni value
pto = nk*pijk[1] + pto; // then convert to pto
@@ -993,7 +993,7 @@ ErrorCode ScdInterface::get_neighbor_sqjk(int np, int pfrom,
facedims[5] = facedims[2];
rdims[5] = ldims[2];
rdims[2] -= dk;
- if (nk < kextra) rdims[2]--;
+ if (pto/pijk[1] < kextra) rdims[2]--;
}
else {
facedims[2] = facedims[5];
@@ -1019,6 +1019,97 @@ ErrorCode ScdInterface::get_neighbor_sqjk(int np, int pfrom,
}
#ifndef USE_MPI
+ErrorCode ScdInterface::get_neighbor_sqijk(int , int ,
+ const int * const , const int * const , const int * const ,
+ int &, int *, int *, int *)
+{
+ return MB_FAILURE;
+#else
+ErrorCode ScdInterface::get_neighbor_sqijk(int np, int pfrom,
+ const int * const gdims, const int * const gperiodic, const int * const dijk,
+ int &pto, int *rdims, int *facedims, int *across_bdy)
+{
+ if (gperiodic[0] || gperiodic[1] || gperiodic[2]) return MB_FAILURE;
+
+ pto = -1;
+ across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
+ int pijk[3], lperiodic[3], ldims[6];
+ ErrorCode rval = compute_partition_sqijk(np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk);
+ if (MB_SUCCESS != rval) return rval;
+ assert(pijk[0] * pijk[1] * pijk[2] == np);
+ pto = -1;
+ bool top[3] = {false, false, false}, bot[3] = {false, false, false};
+ // nijk: rank in i/j/k direction
+ int nijk[3] = {pfrom%pijk[0], (pfrom%(pijk[0]*pijk[1]))/pijk[0], pfrom/(pijk[0]*pijk[1])};
+
+ for (int i = 0; i < 3; i++) {
+ if (nijk[i] == pijk[i]-1) top[i] = true;
+ if (!nijk[i]) bot[i] = true;
+ if ((!gperiodic[i] && bot[i] && -1 == dijk[i]) || // downward && not periodic
+ (!gperiodic[i] && top[i] && 1 == dijk[i])) // upward && not periodic
+ return MB_SUCCESS;
+ }
+
+ std::copy(ldims, ldims+6, facedims);
+ std::copy(ldims, ldims+6, rdims);
+ pto = pfrom;
+ int delijk[3], extra[3];
+ // nijk_to: rank of pto in i/j/k direction
+ int nijk_to[3];
+ for (int i = 0; i < 3; i++) {
+ delijk[i] = (gdims[i+3] == gdims[i] ? 0 : (gdims[i+3] - gdims[i])/pijk[i]);
+ extra[i] = (gdims[i+3]-gdims[i]) % delijk[i];
+ nijk_to[i] = (nijk[i]+dijk[i]+pijk[i]) % pijk[i];
+ }
+ pto = nijk_to[2]*pijk[0]*pijk[1] + nijk_to[1]*pijk[0] + nijk_to[0];
+ assert (pto >= 0 && pto < np);
+ for (int i = 0; i < 3; i++) {
+ if (0 != dijk[i]) {
+ if (-1 == dijk[i]) {
+ facedims[i+3] = facedims[i];
+ if (bot[i]) {
+ // going across lower periodic bdy in i
+ rdims[i+3] = gdims[i+3]+1; // +1 because ldims[4] on remote proc is gdims[4]+1
+ across_bdy[i] = -1;
+ }
+ else {
+ rdims[i+3] = ldims[i];
+ }
+ rdims[i] = rdims[i+3] - delijk[i];
+ if (nijk[i] < extra[i]) rdims[i]--;
+ }
+ else {
+ if (top[i]) {
+ // going across upper periodic bdy in i
+ rdims[i] = gdims[i];
+ facedims[i+3] = gdims[i];
+ across_bdy[i] = 1;
+ }
+ else {
+ rdims[i] = ldims[i+3];
+ }
+ facedims[i] = facedims[i+3];
+ rdims[i+3] = rdims[i] + delijk[i];
+ if (nijk[i] < extra[i]) rdims[i+3]++;
+ if (gperiodic[i] && nijk[i] == dijk[i]-2) rdims[i+3]++; // +1 because next proc is on periodic bdy
+ }
+ }
+ }
+
+ assert(-1 != pto);
+#ifndef NDEBUG
+ for (int i = 0; i < 3; i++) {
+ assert((rdims[i] >= gdims[i] && (rdims[i+3] <= gdims[i+3] || (across_bdy[i] && bot[i]))));
+ assert(((facedims[i] >= rdims[i] || (gperiodic[i] && rdims[i+3] == gdims[i+3] && facedims[i] == gdims[i]))));
+ assert((facedims[i] >= ldims[i] && facedims[i+3] <= ldims[i+3]));
+ }
+#endif
+
+ return MB_SUCCESS;
+#endif
+}
+
+#ifndef USE_MPI
ErrorCode ScdInterface::get_neighbor_alljorkori(int , int ,
const int * const , const int * const , const int * const ,
int &, int *, int *, int *)
@@ -1038,7 +1129,7 @@ ErrorCode ScdInterface::get_neighbor_alljorkori(int np, int pfrom,
if (MB_SUCCESS != rval) return rval;
int ind = -1;
- across_bdy[0] = across_bdy[1] = 0;
+ across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
for (int i = 0; i < 3; i++) {
if (pijk[i] > 1) {
@@ -1126,7 +1217,7 @@ ErrorCode ScdInterface::get_shared_vertices(ParallelComm *pcomm, ScdBox *box,
// get index of partitioned dimension
const int *ldims = box->box_dims();
ErrorCode rval;
- int ijkrem[6], ijkface[6], across_bdy[2];
+ int ijkrem[6], ijkface[6], across_bdy[3];
for (int k = -1; k <= 1; k ++) {
for (int j = -1; j <= 1; j ++) {
diff --git a/src/io/NCHelper.cpp b/src/io/NCHelper.cpp
index b6158a9..c3cd286 100644
--- a/src/io/NCHelper.cpp
+++ b/src/io/NCHelper.cpp
@@ -365,8 +365,8 @@ ErrorCode ScdNCHelper::create_mesh(ScdInterface* scdi, const FileOptions& opts,
std::vector<double>& klVals = _readNC->klVals;
Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
DebugOutput& dbgOut = _readNC->dbgOut;
- int (&locallyPeriodic)[2] = _readNC->locallyPeriodic;
- int (&globallyPeriodic)[2] = _readNC->globallyPeriodic;
+ int (&locallyPeriodic)[3] = _readNC->locallyPeriodic;
+ int (&globallyPeriodic)[3] = _readNC->globallyPeriodic;
ScdParData& parData = _readNC->parData;
Range tmp_range;
diff --git a/src/io/NCHelperEuler.cpp b/src/io/NCHelperEuler.cpp
index 7336fa9..8de0818 100644
--- a/src/io/NCHelperEuler.cpp
+++ b/src/io/NCHelperEuler.cpp
@@ -77,8 +77,8 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
DebugOutput& dbgOut = _readNC->dbgOut;
bool& isParallel = _readNC->isParallel;
int& partMethod = _readNC->partMethod;
- int (&locallyPeriodic)[2] = _readNC->locallyPeriodic;
- int (&globallyPeriodic)[2] = _readNC->globallyPeriodic;
+ int (&locallyPeriodic)[3] = _readNC->locallyPeriodic;
+ int (&globallyPeriodic)[3] = _readNC->globallyPeriodic;
ScdParData& parData = _readNC->parData;
#ifdef USE_MPI
ParallelComm*& myPcomm = _readNC->myPcomm;
@@ -147,7 +147,7 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
#ifdef USE_MPI
for (int i = 0; i < 6; i++)
parData.gDims[i] = gDims[i];
- for (int i = 0; i < 2; i++)
+ for (int i = 0; i < 3; i++)
parData.gPeriodic[i] = globallyPeriodic[i];
parData.partMethod = partMethod;
int pdims[3];
diff --git a/src/io/NCHelperFV.cpp b/src/io/NCHelperFV.cpp
index cb6229e..8b9be3f 100644
--- a/src/io/NCHelperFV.cpp
+++ b/src/io/NCHelperFV.cpp
@@ -72,8 +72,8 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
DebugOutput& dbgOut = _readNC->dbgOut;
bool& isParallel = _readNC->isParallel;
int& partMethod = _readNC->partMethod;
- int (&locallyPeriodic)[2] = _readNC->locallyPeriodic;
- int (&globallyPeriodic)[2] = _readNC->globallyPeriodic;
+ int (&locallyPeriodic)[3] = _readNC->locallyPeriodic;
+ int (&globallyPeriodic)[3] = _readNC->globallyPeriodic;
ScdParData& parData = _readNC->parData;
#ifdef USE_MPI
ParallelComm*& myPcomm = _readNC->myPcomm;
@@ -159,7 +159,7 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
#ifdef USE_MPI
for (int i = 0; i < 6; i++)
parData.gDims[i] = gDims[i];
- for (int i = 0; i < 2; i++)
+ for (int i = 0; i < 3; i++)
parData.gPeriodic[i] = globallyPeriodic[i];
parData.partMethod = partMethod;
int pdims[3];
diff --git a/src/io/ReadNC.hpp b/src/io/ReadNC.hpp
index 882fd98..2909ba2 100644
--- a/src/io/ReadNC.hpp
+++ b/src/io/ReadNC.hpp
@@ -264,10 +264,10 @@ private:
int partMethod;
//! whether mesh is locally periodic in i or j
- int locallyPeriodic[2];
+ int locallyPeriodic[3];
//! whether mesh is globally periodic in i or j
- int globallyPeriodic[2];
+ int globallyPeriodic[3];
//! parallel data object, to be cached with ScdBox
ScdParData parData;
diff --git a/src/io/Tqdcfr.cpp b/src/io/Tqdcfr.cpp
index 833b15b..dc841fc 100644
--- a/src/io/Tqdcfr.cpp
+++ b/src/io/Tqdcfr.cpp
@@ -2491,6 +2491,8 @@ ErrorCode Tqdcfr::parse_acis_attribs(const unsigned int entity_rec_num,
if (records[entity_rec_num].entity == 0) {
records[entity_rec_num].entity = uidSetMap[uid];
}
+
+ assert(records[entity_rec_num].entity);
// set the id
if (id != -1) {
diff --git a/src/moab/ScdInterface.hpp b/src/moab/ScdInterface.hpp
index 3c3e6d1..9612946 100644
--- a/src/moab/ScdInterface.hpp
+++ b/src/moab/ScdInterface.hpp
@@ -106,13 +106,13 @@ class ScdParData {
public:
ScdParData() : partMethod(NOPART) {
gDims[0] = gDims[1] = gDims[2] = gDims[3] = gDims[4] = gDims[5] = -1;
- gPeriodic[0] = gPeriodic[1] = -1;
+ gPeriodic[0] = gPeriodic[1] = gPeriodic[2] = -1;
pDims[0] = pDims[1] = pDims[2] = -1;
}
//! Partition method enumeration; these strategies are described in comments for
//! compute_partition_alljorkori, compute_partition_alljkbal, compute_partition_sqij, and compute_partition_sqjk
- enum PartitionMethod {NOPART=-1, ALLJORKORI=0, ALLJKBAL, SQIJ, SQJK};
+ enum PartitionMethod {NOPART=0, ALLJORKORI, ALLJKBAL, SQIJ, SQJK, SQIJK};
//! partition method used to partition global parametric space
int partMethod;
@@ -121,10 +121,11 @@ public:
int gDims[6];
//! is globally periodic in i or j
- int gPeriodic[2];
+ int gPeriodic[3];
//! number of procs in each direction
int pDims[3];
+
};
class ScdInterface
@@ -325,6 +326,13 @@ private:
inline static ErrorCode compute_partition_sqjk(int np, int nr, const int * const gijk, const int * const gperiodic,
int *lijk, int *lperiodic, int *pijk);
+ //! Compute a partition of structured parameter space
+ /** Partitions the structured parametric space by seeking square ijk partitions
+ * For description of arguments, see ScdInterface::compute_partition.
+ */
+ inline static ErrorCode compute_partition_sqijk(int np, int nr, const int * const gijk, const int * const gperiodic,
+ int *lijk, int *lperiodic, int *pijk);
+
//! Get vertices shared with other processors
/** Shared vertices returned as indices into each proc's handle space
* \param box Box used to get parametric space info
@@ -354,6 +362,10 @@ private:
const int * const gdims, const int * const gperiodic, const int * const dijk,
int &pto, int *rdims, int *facedims, int *across_bdy);
+ static ErrorCode get_neighbor_sqijk(int np, int pfrom,
+ const int * const gdims, const int * const gperiodic, const int * const dijk,
+ int &pto, int *rdims, int *facedims, int *across_bdy);
+
static int gtol(const int *gijk, int i, int j, int k);
//! interface instance
@@ -597,11 +609,17 @@ public:
*/
bool locally_periodic_j() const;
+ //! Return whether box is locally periodic in k
+ /** Return whether box is locally periodic in k
+ * \return True if box is locally periodic in k direction
+ */
+ bool locally_periodic_k() const;
+
//! Return whether box is locally periodic in i and j
/** Return whether box is locally periodic in i and j
* \param lperiodic Non-zero if locally periodic in i [0] or j [1]
*/
- void locally_periodic(bool lperiodic[2]) const;
+ void locally_periodic(bool lperiodic[3]) const;
//! Return parallel data
/** Return parallel data, if there is any
@@ -685,7 +703,7 @@ private:
int boxDims[6];
//! is locally periodic in i or j
- int locallyPeriodic[2];
+ int locallyPeriodic[3];
//! parallel data associated with this box, if any
ScdParData parData;
@@ -722,6 +740,10 @@ inline ErrorCode ScdInterface::compute_partition(int np, int nr, const ScdParDat
rval = compute_partition_sqjk(np, nr, par_data.gDims, par_data.gPeriodic,
ldims, lperiodic, pdims);
break;
+ case ScdParData::SQIJK:
+ rval = compute_partition_sqijk(np, nr, par_data.gDims, par_data.gPeriodic,
+ ldims, lperiodic, pdims);
+ break;
default:
rval = MB_FAILURE;
break;
@@ -740,8 +762,7 @@ inline ErrorCode ScdInterface::compute_partition_alljorkori(int np, int nr,
if (!lperiodic) lperiodic = tmp_lp;
if (!pijk) pijk = tmp_pijk;
- lperiodic[0] = gperiodic[0];
- lperiodic[1] = gperiodic[1];
+ for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];
if (np == 1) {
if (ldims) {
@@ -822,8 +843,7 @@ inline ErrorCode ScdInterface::compute_partition_alljkbal(int np, int nr,
if (!lperiodic) lperiodic = tmp_lp;
if (!pijk) pijk = tmp_pijk;
- lperiodic[0] = gperiodic[0];
- lperiodic[1] = gperiodic[1];
+ for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];
if (np == 1) {
if (ldims) {
@@ -899,8 +919,7 @@ inline ErrorCode ScdInterface::compute_partition_sqij(int np, int nr,
// square IxJ partition
- lperiodic[0] = gperiodic[0];
- lperiodic[1] = gperiodic[1];
+ for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];
if (np == 1) {
if (ldims) {
@@ -979,8 +998,7 @@ inline ErrorCode ScdInterface::compute_partition_sqjk(int np, int nr,
if (!pijk) pijk = tmp_pijk;
// square JxK partition
- lperiodic[0] = gperiodic[0];
- lperiodic[1] = gperiodic[1];
+ for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];
if (np == 1) {
if (ldims) {
@@ -1042,6 +1060,91 @@ inline ErrorCode ScdInterface::compute_partition_sqjk(int np, int nr,
return MB_SUCCESS;
}
+inline ErrorCode ScdInterface::compute_partition_sqijk(int np, int nr,
+ const int * const gijk, const int * const gperiodic,
+ int *ldims, int *lperiodic, int *pijk)
+{
+ if (gperiodic[0] || gperiodic[1] || gperiodic[2]) return MB_FAILURE;
+
+ int tmp_lp[3], tmp_pijk[3];
+ if (!lperiodic) lperiodic = tmp_lp;
+ if (!pijk) pijk = tmp_pijk;
+
+ // square IxJxK partition
+
+ for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];
+
+ if (np == 1) {
+ if (ldims)
+ for (int i = 0; i < 6; i++) ldims[i] = gijk[i];
+ pijk[0] = pijk[1] = pijk[2] = 1;
+ return MB_SUCCESS;
+ }
+
+ std::vector<int> pfactors;
+ pfactors.push_back(1);
+ for (int i = 2; i <= np/2; i++)
+ if (!(np%i)) pfactors.push_back(i);
+ pfactors.push_back(np);
+
+ // test for IJ, JK, IK
+ int IJK[3], dIJK[3];
+ for (int i = 0; i < 3; i++)
+ IJK[i] = std::max(gijk[3+i] - gijk[i], 1);
+ // order IJK from lo to hi
+ int lo = 0, hi = 0;
+ for (int i = 1; i < 3; i++) {
+ if (IJK[i] < IJK[lo]) lo = i;
+ if (IJK[i] > IJK[hi]) hi = i;
+ }
+ if (lo == hi) hi = (lo+1)%3;
+ int mid = 3 - lo - hi;
+ // search for perfect subdivision of np that balances #cells
+ int perfa_best = -1, perfb_best = -1;
+ double ratio = 0.0;
+ for (int po = 0; po < (int)pfactors.size(); po++) {
+ for (int pi = po; pi < (int)pfactors.size() && np/(pfactors[po]*pfactors[pi]) >= pfactors[pi]; pi++) {
+ int p3 = std::find(pfactors.begin(), pfactors.end(), np/(pfactors[po]*pfactors[pi])) - pfactors.begin();
+ if (p3 == (int)pfactors.size() || pfactors[po]*pfactors[pi]*pfactors[p3] != np) continue; // po*pi should exactly factor np
+ assert(po <= pi && pi <= p3);
+ // by definition, po <= pi <= p3
+ double minl = std::min(std::min(IJK[lo]/pfactors[po], IJK[mid]/pfactors[pi]), IJK[hi]/pfactors[p3]),
+ maxl = std::max(std::max(IJK[lo]/pfactors[po], IJK[mid]/pfactors[pi]), IJK[hi]/pfactors[p3]);
+ if (minl/maxl > ratio) {
+ ratio = minl/maxl;
+ perfa_best = po; perfb_best = pi;
+ }
+ }
+ }
+ if (perfa_best == -1 || perfb_best == -1) return MB_FAILURE;
+
+ // VARIABLES DESCRIBING THE MESH:
+ // pijk[i] = # procs in direction i
+ // numr[i] = my proc's position in direction i
+ // dIJK[i] = # edges/elements in direction i
+ // extra[i]= # procs having extra edge in direction i
+ // top[i] = if true, I'm the last proc in direction i
+
+ pijk[lo] = pfactors[perfa_best]; pijk[mid] = pfactors[perfb_best];
+ pijk[hi] = (np/(pfactors[perfa_best]*pfactors[perfb_best]));
+ int extra[3] = {0, 0, 0}, numr[3];
+ for (int i = 0; i < 3; i++) {
+ dIJK[i] = IJK[i] / pijk[i];
+ extra[i] = IJK[i]%pijk[i];
+ }
+ numr[2] = nr / (pijk[0]*pijk[1]);
+ int rem = nr % (pijk[0] * pijk[1]);
+ numr[1] = rem / pijk[0];
+ numr[0] = rem % pijk[0];
+ for (int i = 0; i < 3; i++) {
+ extra[i] = IJK[i] % dIJK[i];
+ ldims[i] = gijk[i] + numr[i]*dIJK[i] + std::min(extra[i], numr[i]);
+ ldims[3+i] = ldims[i] + dIJK[i] + (numr[i] < extra[i] ? 1 : 0);
+ }
+
+ return MB_SUCCESS;
+}
+
inline int ScdInterface::gtol(const int *gijk, int i, int j, int k)
{
return ((k-gijk[2])*(gijk[3]-gijk[0]+1)*(gijk[4]-gijk[1]+1) + (j-gijk[1])*(gijk[3]-gijk[0]+1) + i-gijk[0]);
@@ -1097,6 +1200,9 @@ inline ErrorCode ScdInterface::get_neighbor(int np, int pfrom, const ScdParData
case ScdParData::SQJK:
return get_neighbor_sqjk(np, pfrom, spd.gDims, spd.gPeriodic, dijk,
pto, rdims, facedims, across_bdy);
+ case ScdParData::SQIJK:
+ return get_neighbor_sqijk(np, pfrom, spd.gDims, spd.gPeriodic, dijk,
+ pto, rdims, facedims, across_bdy);
default:
break;
}
@@ -1134,7 +1240,7 @@ inline int ScdBox::num_elements() const
return (!startElem ? 0 :
(boxSize[0]- (locallyPeriodic[0] ? 0 : 1)) *
(-1 == boxSize[1] ? 1 : (boxSize[1]-(locallyPeriodic[1] ? 0 : 1))) *
- (boxSize[2] == -1 || boxSize[2] == 1 ? 1 : (boxSize[2]-1)));
+ (boxSize[2] == -1 || boxSize[2] == 1 ? 1 : (boxSize[2]-(locallyPeriodic[2] ? 0 : 1))));
}
inline int ScdBox::num_vertices() const
@@ -1267,11 +1373,27 @@ inline bool ScdBox::locally_periodic_j() const
return locallyPeriodic[1];
}
-inline void ScdBox::locally_periodic(bool lperiodic[2]) const
+inline bool ScdBox::locally_periodic_k() const
+{
+ return locallyPeriodic[2];
+}
+
+inline void ScdBox::locally_periodic(bool lperiodic[3]) const
{
- for (int i = 0; i < 2; i++)
+ for (int i = 0; i < 3; i++)
lperiodic[i] = locallyPeriodic[i];
}
+inline std::ostream &operator<<(std::ostream &str, const ScdParData &pd)
+{
+ static const char *PartitionMethodNames[] = {"NOPART", "ALLJORKORI", "ALLJKBAL", "SQIJ", "SQJK", "SQIJK"};
+ str << "Partition method = " << PartitionMethodNames[pd.partMethod] << ", gDims = ("
+ << pd.gDims[0] << "," << pd.gDims[1] << "," << pd.gDims[2] << ")-("
+ << pd.gDims[3] << "," << pd.gDims[4] << "," << pd.gDims[5] << "), gPeriodic = ("
+ << pd.gPeriodic[0] << ", " << pd.gPeriodic[1] << "," << pd.gPeriodic[2] << "), pDims = ("
+ << pd.pDims[0] << "," << pd.pDims[1] << "," << pd.pDims[2] << ")" << std::endl;
+ return str;
+}
+
} // namespace moab
#endif
diff --git a/test/Makefile.am b/test/Makefile.am
index fbf76e1..5b63e9e 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -27,6 +27,7 @@ TESTS = range_test \
moab_test \
homxform_test \
scdseq_test \
+ scd_test_partn \
test_adj \
seq_man_test \
geom_util_test \
@@ -87,6 +88,7 @@ homxform_test_SOURCES = $(top_srcdir)/src/HomXform.cpp
homxform_test_CPPFLAGS = -DTEST $(AM_CPPFLAGS) $(CPPFLAGS)
scdseq_test_SOURCES = scdseq_test.cpp
+scd_test_partn_SOURCES = scd_test_partn.cpp
#scdseq_timing_SOURCES = scdseq_timing.cpp
#merge_test_SOURCES = merge_test.cpp
#test_exo_SOURCES = test_exo.cpp
diff --git a/test/scd_test_partn.cpp b/test/scd_test_partn.cpp
new file mode 100644
index 0000000..59a63c5
--- /dev/null
+++ b/test/scd_test_partn.cpp
@@ -0,0 +1,68 @@
+#include "moab/ScdInterface.hpp"
+#include "moab/Core.hpp"
+#include "TestUtil.hpp"
+
+#include <iostream>
+
+#ifdef USE_MPI
+#include "moab/ParallelComm.hpp"
+#endif
+
+using namespace moab;
+
+int gdims[6], np;
+
+void test_sqijk()
+{
+ Core moab;
+ ScdInterface *scd;
+ ErrorCode rval = moab.Interface::query_interface(scd);
+ CHECK_ERR(rval);
+ int ldims[6];
+ ScdParData par_data;
+ std::copy(gdims, gdims+6, par_data.gDims);
+ for (int i = 0; i < 3; i++) par_data.gPeriodic[i] = 0;
+ par_data.partMethod = ScdParData::SQIJK;
+
+ std::cout << "gDims = (" << par_data.gDims[0] << "," << par_data.gDims[1] << "," << par_data.gDims[2] << ")--("
+ << par_data.gDims[3] << "," << par_data.gDims[4] << "," << par_data.gDims[5] << ")" << std::endl;
+
+ int lperiodic[3], pijk[3];
+
+ rval = ScdInterface::compute_partition(np, 0, par_data, ldims, lperiodic, pijk); CHECK_ERR(rval);
+
+
+ std::cout << "#proc in 3 directions = (" << pijk[0] << "," << pijk[1] << "," << pijk[2] << ")" << std::endl;
+ std::cout << "local dims are ("
+ << ldims[0] << "," << ldims[1] << "," << ldims[2] << ")--("
+ << ldims[3] << "," << ldims[4] << "," << ldims[5] << ")\n";
+}
+
+int main(int argc, char**argv)
+{
+ if (argc < 2) {
+ std::cout << "Usage: " << argv[0] << " <#proc> [<imax> [<jmax><kmax>]]" << std::endl;
+ exit(1);
+ }
+
+ if (argc < 3) {
+ np = atoi(argv[1]);
+ gdims[0] = gdims[1] = gdims[2] = 0;
+ gdims[3] = gdims[4] = gdims[5] = 100;
+ }
+ else if (argc < 4) {
+ np = atoi(argv[1]);
+ gdims[0] = gdims[1] = gdims[2] = 0;
+ gdims[3] = gdims[4] = gdims[5] = atoi(argv[2]);
+ }
+ else if (argc < 6) {
+ np = atoi(argv[1]);
+ gdims[0] = gdims[1] = gdims[2] = 0;
+ gdims[3] = atoi(argv[2]);
+ gdims[4] = atoi(argv[3]);
+ gdims[5] = atoi(argv[4]);
+ }
+
+ // test partition method
+ RUN_TEST(test_sqijk);
+}
diff --git a/test/scdseq_test.cpp b/test/scdseq_test.cpp
index 651eec9..e5acb0d 100644
--- a/test/scdseq_test.cpp
+++ b/test/scdseq_test.cpp
@@ -264,17 +264,18 @@ ErrorCode evaluate_element_sequence(ScdBox *this_box)
int main(int, char**)
{
// test partition methods
- RUN_TEST(test_parallel_partitions);
+ int err = RUN_TEST(test_parallel_partitions);
// test creating and evaluating vertex sequences
- RUN_TEST(test_vertex_seq);
+ err += RUN_TEST(test_vertex_seq);
// test creating and evaluating element sequences
- RUN_TEST(test_element_seq);
+ err += RUN_TEST(test_element_seq);
// test periodic sequences
- RUN_TEST(test_periodic_seq);
+ err += RUN_TEST(test_periodic_seq);
+ return err;
}
void test_vertex_seq()
@@ -1324,38 +1325,48 @@ void test_parallel_partitions()
#ifndef NDEBUG
maxpow = 4;
#endif
+
+ int fails = 0;
for (int exp = 2; exp <= maxpow; exp += 2) {
int nprocs = 0.1 + pow(2.0, (double)exp);
// alljorkori
rval = test_parallel_partition(gdims, nprocs, ScdParData::ALLJORKORI);
- CHECK_ERR(rval);
+ if (MB_SUCCESS != rval) fails++;
// alljkbal
rval = test_parallel_partition(gdims, nprocs, ScdParData::ALLJKBAL);
- CHECK_ERR(rval);
+ if (MB_SUCCESS != rval) fails++;
// sqij
rval = test_parallel_partition(gdims, nprocs, ScdParData::SQIJ);
- CHECK_ERR(rval);
+ if (MB_SUCCESS != rval) fails++;
// sqjk
rval = test_parallel_partition(gdims, nprocs, ScdParData::SQJK);
- CHECK_ERR(rval);
+ if (MB_SUCCESS != rval) fails++;
+
+ // sqijk
+ rval = test_parallel_partition(gdims, nprocs, ScdParData::SQIJK);
+ if (MB_SUCCESS != rval) fails++;
}
+
+ if (fails) CHECK_ERR(MB_FAILURE);
}
ErrorCode test_parallel_partition(int *gdims, int nprocs, int part_method)
{
ErrorCode rval;
- int pto, pfrom, across_bdy_a[2], across_bdy_b[2], rdims_a[6], rdims_b[6], facedims_a[6], facedims_b[6], ldims[6];
+ int pto, pfrom, across_bdy_a[2], across_bdy_b[2], rdims_a[6], rdims_b[6], facedims_a[6], facedims_b[6], ldims[6],
+ lper[2], pijk[3];
ScdParData spd;
for (int i = 0; i < 6; i++) spd.gDims[i] = gdims[i];
- for (int i = 0; i < 2; i++) spd.gPeriodic[i] = 0;
+ for (int i = 0; i < 3; i++) spd.gPeriodic[i] = 0;
spd.partMethod = part_method;
+ int fails = 0;
for (int p = 0; p < nprocs/2; p++) {
- rval = ScdInterface::compute_partition(nprocs, p, spd, ldims);
+ rval = ScdInterface::compute_partition(nprocs, p, spd, ldims, lper, pijk);
if (MB_SUCCESS != rval) continue;
for (int k = -1; k <= 1; k++) {
@@ -1367,9 +1378,10 @@ ErrorCode test_parallel_partition(int *gdims, int nprocs, int part_method)
if (MB_SUCCESS != rval) return rval;
if (-1 == pto) continue;
+ bool fail = false;
if (facedims_a[0] < rdims_a[0] || facedims_a[0] > rdims_a[3] ||
facedims_a[1] < rdims_a[1] || facedims_a[1] > rdims_a[4] ||
- facedims_a[2] < rdims_a[2] || facedims_a[2] > rdims_a[5]) CHECK_ERR(MB_FAILURE);
+ facedims_a[2] < rdims_a[2] || facedims_a[2] > rdims_a[5]) fail = true;
// non-negative value of pto; check corresponding input from that proc to this
@@ -1378,18 +1390,29 @@ ErrorCode test_parallel_partition(int *gdims, int nprocs, int part_method)
pfrom, rdims_b, facedims_b, across_bdy_b);
if (MB_SUCCESS != rval) return rval;
for (int ind = 0; ind < 3; ind++)
- if (facedims_a[ind] < rdims_b[ind] || facedims_b[ind] > rdims_b[ind+3]) CHECK_ERR(MB_FAILURE);
+ if (facedims_a[ind] < rdims_b[ind] || facedims_b[ind] > rdims_b[ind+3]) fail = true;
for (int ind = 0; ind < 6; ind++) {
- if (facedims_a[ind] != facedims_b[ind]) CHECK_ERR(MB_FAILURE);
- if (rdims_b[ind] != ldims[ind]) CHECK_ERR(MB_FAILURE);
+ if (facedims_a[ind] != facedims_b[ind] || rdims_b[ind] != ldims[ind]) fail = true;
+ }
+ if (across_bdy_a[0] != across_bdy_b[0] || across_bdy_a[1] != across_bdy_b[1]) fail = true;
+#define PARRAY(a) "(" << a[0] << "," << a[1] << "," << a[2] << ")"
+#define PARRAY3(a,b,c) "(" << a << "," << b << "," << c << ")"
+#define PARRAY6(a) PARRAY(a) << "-" << PARRAY((a+3))
+ if (fail) {
+ fails++;
+ for (int l = 0; l < 3; l++) spd.pDims[l] = pijk[l];
+ std::cerr << "partMeth = " << part_method << ", p/np = " << p << "/" << nprocs
+ << ", (i,j,k) = " << PARRAY3(i,j,k) << ", ldims = " << PARRAY6(ldims) << std::endl
+ << "facedims_a = " << PARRAY6(facedims_a)
+ << ", facedims_b = " << PARRAY6(facedims_b) << std::endl
+ << "rdims_a = " << PARRAY6(rdims_a) << ", rdims_b = " << PARRAY6(rdims_b) << std::endl;
+ std::cerr << "ScdParData: " << spd << std::endl;
}
-
- if (across_bdy_a[0] != across_bdy_b[0] || across_bdy_a[1] != across_bdy_b[1]) CHECK_ERR(MB_FAILURE);
} // i
} // j
} // k
} // p
- return MB_SUCCESS;
+ return (fails ? MB_FAILURE : MB_SUCCESS);
}
https://bitbucket.org/fathomteam/moab/commits/0a8790877bff/
Changeset: 0a8790877bff
Branch: master
User: tautges
Date: 2013-08-12 21:53:17
Summary: Merge branch 'master' of git at bitbucket.org:fathomteam/moab.git
Conflicts:
src/io/NCHelper.cpp
src/io/NCHelperEuler.cpp
src/io/NCHelperFV.cpp
src/io/ReadNC.hpp
Affected #: 9 files
diff --git a/src/io/NCHelper.cpp b/src/io/NCHelper.cpp
index c3cd286..851d9b1 100644
--- a/src/io/NCHelper.cpp
+++ b/src/io/NCHelper.cpp
@@ -4,6 +4,8 @@
#include "NCHelperHOMME.hpp"
#include "NCHelperMPAS.hpp"
+#include <sstream>
+
#include "moab/ReadUtilIface.hpp"
#include "MBTagConventions.hpp"
@@ -52,89 +54,254 @@ NCHelper* NCHelper::get_nc_helper(ReadNC* readNC, int fileId, const FileOptions&
return NULL;
}
-ErrorCode NCHelper::read_variable_to_set_allocate(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums)
-{
+ErrorCode NCHelper::create_conventional_tags(ScdInterface* scdi, EntityHandle file_set, const std::vector<int>& tstep_nums) {
+ Interface*& mbImpl = _readNC->mbImpl;
+ std::vector<std::string>& dimNames = _readNC->dimNames;
std::vector<int>& dimVals = _readNC->dimVals;
- int tDim = _readNC->tDim;
+ std::map<std::string, ReadNC::AttData>& globalAtts = _readNC->globalAtts;
+ std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
DebugOutput& dbgOut = _readNC->dbgOut;
+ int& partMethod = _readNC->partMethod;
+
+ ErrorCode rval;
+ std::string tag_name;
+
+ // <__NUM_DIMS>
+ Tag numDimsTag = 0;
+ tag_name = "__NUM_DIMS";
+ int numDims = dimNames.size();
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 1, MB_TYPE_INTEGER, numDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __NUM_DIMS tag.");
+ rval = mbImpl->tag_set_data(numDimsTag, &file_set, 1, &numDims);
+ ERRORR(rval, "Trouble setting data for __NUM_DIMS tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+
+ // <__NUM_VARS>
+ Tag numVarsTag = 0;
+ tag_name = "__NUM_VARS";
+ int numVars = varInfo.size();
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 1, MB_TYPE_INTEGER, numVarsTag, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __NUM_VARS tag.");
+ rval = mbImpl->tag_set_data(numVarsTag, &file_set, 1, &numVars);
+ ERRORR(rval, "Trouble setting data for __NUM_VARS tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+
+ // <__DIM_NAMES>
+ Tag dimNamesTag = 0;
+ tag_name = "__DIM_NAMES";
+ std::string dimnames;
+ unsigned int dimNamesSz = dimNames.size();
+ for (unsigned int i = 0; i != dimNamesSz; ++i) {
+ dimnames.append(dimNames[i]);
+ dimnames.push_back('\0');
+ }
+ int dimnamesSz = dimnames.size();
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, dimNamesTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating __DIM_NAMES tag.");
+ const void* ptr = dimnames.c_str();
+ rval = mbImpl->tag_set_by_ptr(dimNamesTag, &file_set, 1, &ptr, &dimnamesSz);
+ ERRORR(rval, "Trouble setting data for __DIM_NAMES tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+
+ // <__DIM_VALUES>
+ Tag dimValsTag = 0;
+ tag_name = "__DIM_VALUES";
+ int dimValsSz = (int)dimVals.size();
+
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_INTEGER, dimValsTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating __DIM_VALUES tag.");
+ ptr = &(dimVals[0]);
+ rval = mbImpl->tag_set_by_ptr(dimValsTag, &file_set, 1, &ptr, &dimValsSz);
+ ERRORR(rval, "Trouble setting data for __DIM_VALUES tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+
+ // <__VAR_NAMES>
+ Tag varNamesTag = 0;
+ tag_name = "__VAR_NAMES";
+ std::string varnames;
+ std::map<std::string, ReadNC::VarData>::iterator mapIter;
+ for (mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter) {
+ varnames.append(mapIter->first);
+ varnames.push_back('\0');
+ }
+ int varnamesSz = varnames.size();
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, varNamesTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating __VAR_NAMES tag.");
+ ptr = varnames.c_str();
+ rval = mbImpl->tag_set_by_ptr(varNamesTag, &file_set, 1, &ptr, &varnamesSz);
+ ERRORR(rval, "Trouble setting data for __VAR_NAMES tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+
+ // __<dim_name>_LOC_MINMAX
+ for (unsigned int i = 0; i != dimNamesSz; ++i) {
+ if (dimNames[i] == "time") {
+ std::stringstream ss_tag_name;
+ ss_tag_name << "__" << dimNames[i] << "_LOC_MINMAX";
+ tag_name = ss_tag_name.str();
+ Tag tagh = 0;
+ std::vector<int> val(2, 0);
+ val[0] = 0;
+ val[1] = nTimeSteps - 1;
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<dim_name>_LOC_MINMAX tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<dim_name>_LOC_MINMAX tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+ }
- ErrorCode rval = MB_SUCCESS;
-
- for (unsigned int i = 0; i < vdatas.size(); i++) {
- if ((std::find(vdatas[i].varDims.begin(), vdatas[i].varDims.end(), tDim) != vdatas[i].varDims.end()))
- vdatas[i].has_t = true;
-
- for (unsigned int t = 0; t < tstep_nums.size(); t++) {
- dbgOut.tprintf(2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
-
- // Get the tag to read into
- if (!vdatas[i].varTags[t]) {
- rval = _readNC->get_tag_to_set(vdatas[i], tstep_nums[t], vdatas[i].varTags[t]);
- ERRORR(rval, "Trouble getting tag.");
- }
+ // __<dim_name>_LOC_VALS
+ for (unsigned int i = 0; i != dimNamesSz; ++i) {
+ if (dimNames[i] != "time")
+ continue;
+ std::vector<int> val;
+ if (!tstep_nums.empty())
+ val = tstep_nums;
+ else {
+ val.resize(tVals.size());
+ for (unsigned int j = 0; j != tVals.size(); ++j)
+ val[j] = j;
+ }
+ Tag tagh = 0;
+ std::stringstream ss_tag_name;
+ ss_tag_name << "__" << dimNames[i] << "_LOC_VALS";
+ tag_name = ss_tag_name.str();
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), val.size(), MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<dim_name>_LOC_VALS tag.");
+ rval = mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<dim_name>_LOC_VALS tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
- // Assume point-based values for now?
- if (-1 == tDim || dimVals[tDim] <= (int) t)
- ERRORR(MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number.");
+ // __<var_name>_DIMS
+ for (mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter) {
+ Tag varNamesDimsTag = 0;
+ std::stringstream ss_tag_name;
+ ss_tag_name << "__" << mapIter->first << "_DIMS";
+ tag_name = ss_tag_name.str();
+ unsigned int varDimSz = varInfo[mapIter->first].varDims.size();
+ if (varDimSz == 0)
+ continue;
+ varInfo[mapIter->first].varTags.resize(varDimSz, 0);
+ for (unsigned int i = 0; i != varDimSz; ++i) {
+ Tag tmptag = 0;
+ std::string tmptagname = dimNames[varInfo[mapIter->first].varDims[i]];
+ mbImpl->tag_get_handle(tmptagname.c_str(), 0, MB_TYPE_OPAQUE, tmptag, MB_TAG_ANY);
+ varInfo[mapIter->first].varTags[i] = tmptag;
+ }
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), varDimSz, MB_TYPE_HANDLE, varNamesDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<var_name>_DIMS tag.");
+ rval = mbImpl->tag_set_data(varNamesDimsTag, &file_set, 1, &(varInfo[mapIter->first].varTags[0]));
+ ERRORR(rval, "Trouble setting data for __<var_name>_DIMS tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
- // Set up the dimensions and counts
- // First variable dimension is time, if it exists
- if (vdatas[i].has_t)
- {
- if (vdatas[i].varDims.size() != 1)
- {
- vdatas[i].readStarts[t].push_back(tstep_nums[t]);
- vdatas[i].readCounts[t].push_back(1);
- }
- else
- {
- vdatas[i].readStarts[t].push_back(0);
- vdatas[i].readCounts[t].push_back(tstep_nums.size());
- }
- }
+ // <PARTITION_METHOD>
+ Tag part_tag = scdi->part_method_tag();
+ if (!part_tag)
+ ERRORR(MB_FAILURE, "Trouble getting partition method tag.");
+ rval = mbImpl->tag_set_data(part_tag, &file_set, 1, &partMethod);
+ ERRORR(rval, "Trouble setting data for PARTITION_METHOD tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+
+ // <__GLOBAL_ATTRIBS>
+ tag_name = "__GLOBAL_ATTRIBS";
+ Tag globalAttTag = 0;
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, globalAttTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating __GLOBAL_ATTRIBS tag.");
+ std::string gattVal;
+ std::vector<int> gattLen;
+ rval = create_attrib_string(globalAtts, gattVal, gattLen);
+ ERRORR(rval, "Trouble creating attribute strings.");
+ const void* gattptr = gattVal.c_str();
+ int globalAttSz = gattVal.size();
+ rval = mbImpl->tag_set_by_ptr(globalAttTag, &file_set, 1, &gattptr, &globalAttSz);
+ ERRORR(rval, "Trouble setting data for __GLOBAL_ATTRIBS tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+
+ // <__GLOBAL_ATTRIBS_LEN>
+ tag_name = "__GLOBAL_ATTRIBS_LEN";
+ Tag globalAttLenTag = 0;
+ if (gattLen.size() == 0)
+ gattLen.push_back(0);
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), gattLen.size(), MB_TYPE_INTEGER, globalAttLenTag, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __GLOBAL_ATTRIBS_LEN tag.");
+ rval = mbImpl->tag_set_data(globalAttLenTag, &file_set, 1, &gattLen[0]);
+ ERRORR(rval, "Trouble setting data for __GLOBAL_ATTRIBS_LEN tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+
+ // __<var_name>_ATTRIBS and __<var_name>_ATTRIBS_LEN
+ for (mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter) {
+ std::stringstream ssTagName;
+ ssTagName << "__" << mapIter->first << "_ATTRIBS";
+ tag_name = ssTagName.str();
+ Tag varAttTag = 0;
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, varAttTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating __<var_name>_ATTRIBS tag.");
+ std::string varAttVal;
+ std::vector<int> varAttLen;
+ rval = create_attrib_string(mapIter->second.varAtts, varAttVal, varAttLen);
+ ERRORR(rval, "Trouble creating attribute strings.");
+ const void* varAttPtr = varAttVal.c_str();
+ int varAttSz = varAttVal.size();
+ rval = mbImpl->tag_set_by_ptr(varAttTag, &file_set, 1, &varAttPtr, &varAttSz);
+ ERRORR(rval, "Trouble setting data for __<var_name>_ATTRIBS tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ if (varAttLen.size() == 0)
+ varAttLen.push_back(0);
+ ssTagName << "_LEN";
+ tag_name = ssTagName.str();
+ Tag varAttLenTag = 0;
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), varAttLen.size(), MB_TYPE_INTEGER, varAttLenTag, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<var_name>_ATTRIBS_LEN tag.");
+ rval = mbImpl->tag_set_data(varAttLenTag, &file_set, 1, &varAttLen[0]);
+ ERRORR(rval, "Trouble setting data for __<var_name>_ATTRIBS_LEN tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
- // Set up other dimensions and counts
- if (vdatas[i].varDims.empty()) {
- // Scalar variable
- vdatas[i].readStarts[t].push_back(0);
- vdatas[i].readCounts[t].push_back(1);
- }
- else {
- for (unsigned int idx = 0; idx != vdatas[i].varDims.size(); idx++){
- if (tDim != vdatas[i].varDims[idx]){
- // Push other variable dimensions, except time, which was already pushed
- vdatas[i].readStarts[t].push_back(0);
- vdatas[i].readCounts[t].push_back(dimVals[vdatas[i].varDims[idx]]);
- }
- }
- }
- std::size_t sz = 1;
- for (std::size_t idx = 0; idx != vdatas[i].readCounts[t].size(); idx++)
- sz *= vdatas[i].readCounts[t][idx];
- vdatas[i].sz = sz;
- switch (vdatas[i].varDataType) {
- case NC_BYTE:
- case NC_CHAR:
- vdatas[i].varDatas[t] = new char[sz];
- break;
- case NC_DOUBLE:
- case NC_FLOAT:
- vdatas[i].varDatas[t] = new double[sz];
- break;
- case NC_INT:
- case NC_SHORT:
- vdatas[i].varDatas[t] = new int[sz];
- break;
- default:
- std::cerr << "Unrecognized data type for tag " << std::endl;
- rval = MB_FAILURE;
- }
- if (vdatas[i].varDims.size() <= 1)
- break;
- }
+ // <__VAR_NAMES_LOCATIONS>
+ tag_name = "__VAR_NAMES_LOCATIONS";
+ Tag varNamesLocsTag = 0;
+ std::vector<int> varNamesLocs(varInfo.size());
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), varNamesLocs.size(), MB_TYPE_INTEGER, varNamesLocsTag, MB_TAG_CREAT
+ | MB_TAG_SPARSE);
+ ERRORR(rval, "Trouble creating __VAR_NAMES_LOCATIONS tag.");
+ for (mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter) {
+ varNamesLocs[std::distance(varInfo.begin(), mapIter)] = mapIter->second.entLoc;
}
+ rval = mbImpl->tag_set_data(varNamesLocsTag, &file_set, 1, &varNamesLocs[0]);
+ ERRORR(rval, "Trouble setting data for __VAR_NAMES_LOCATIONS tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+
+ // <__MESH_TYPE>
+ Tag meshTypeTag = 0;
+ tag_name = "__MESH_TYPE";
+ std::string meshTypeName = 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.");
+ ptr = meshTypeName.c_str();
+ int leng= meshTypeName.size();
+ rval = mbImpl->tag_set_by_ptr(meshTypeTag, &file_set, 1, &ptr, &leng);
+ ERRORR(rval, "Trouble setting data for __MESH_TYPE tag.");
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
- return rval;
+ return MB_SUCCESS;
}
ErrorCode NCHelper::read_variable_to_set(EntityHandle file_set, std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums)
@@ -314,9 +481,301 @@ ErrorCode NCHelper::convert_variable(ReadNC::VarData& var_data, int tstep_num)
return MB_SUCCESS;
}
+ErrorCode NCHelper::read_coordinate(const char* var_name, int lmin, int lmax, std::vector<double>& cvals)
+{
+ std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
+ std::map<std::string, ReadNC::VarData>::iterator vmit = varInfo.find(var_name);
+ if (varInfo.end() == vmit)
+ return MB_FAILURE;
+
+ assert(lmin >= 0 && lmax >= lmin);
+ NCDF_SIZE tstart = lmin;
+ NCDF_SIZE tcount = lmax - lmin + 1;
+ NCDF_DIFF dum_stride = 1;
+ int fail;
+
+ // Check size
+ if (tcount != cvals.size())
+ cvals.resize(tcount);
+
+ // Check to make sure it's a float or double
+ if (NC_DOUBLE == (*vmit).second.varDataType) {
+ fail = NCFUNCA(get_vars_double)(_fileId, (*vmit).second.varId, &tstart, &tcount, &dum_stride, &cvals[0]);
+ if (fail)
+ ERRORS(MB_FAILURE, "Failed to get coordinate values.");
+ }
+ else if (NC_FLOAT == (*vmit).second.varDataType) {
+ std::vector<float> tcvals(tcount);
+ fail = NCFUNCA(get_vars_float)(_fileId, (*vmit).second.varId, &tstart, &tcount, &dum_stride, &tcvals[0]);
+ if (fail)
+ ERRORS(MB_FAILURE, "Failed to get coordinate values.");
+ std::copy(tcvals.begin(), tcvals.end(), cvals.begin());
+ }
+ else {
+ ERRORR(MB_FAILURE, "Wrong data type for coordinate variable.");
+ }
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelper::get_tag_to_set(ReadNC::VarData& var_data, int tstep_num, Tag& tagh)
+{
+ Interface*& mbImpl = _readNC->mbImpl;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+
+ std::ostringstream tag_name;
+ if ((!var_data.has_t) || (var_data.varDims.size() <= 1))
+ tag_name << var_data.varName;
+ else if (!tstep_num) {
+ std::string tmp_name = var_data.varName + "0";
+ tag_name << tmp_name.c_str();
+ }
+ else
+ tag_name << var_data.varName << tstep_num;
+ ErrorCode rval = MB_SUCCESS;
+ tagh = 0;
+ switch (var_data.varDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ rval = mbImpl->tag_get_handle(tag_name.str().c_str(), 0, MB_TYPE_OPAQUE, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ break;
+ case NC_DOUBLE:
+ case NC_FLOAT:
+ rval = mbImpl->tag_get_handle(tag_name.str().c_str(), 0, MB_TYPE_DOUBLE, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ break;
+ case NC_INT:
+ case NC_SHORT:
+ rval = mbImpl->tag_get_handle(tag_name.str().c_str(), 0, MB_TYPE_INTEGER, tagh, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ break;
+ default:
+ std::cerr << "Unrecognized data type for tag " << tag_name << std::endl;
+ rval = MB_FAILURE;
+ }
+
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.str().c_str());
+
+ return rval;
+}
+
+ErrorCode NCHelper::get_tag_to_nonset(ReadNC::VarData& var_data, int tstep_num, Tag& tagh, int num_lev)
+{
+ Interface*& mbImpl = _readNC->mbImpl;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+
+ std::ostringstream tag_name;
+ if (!tstep_num) {
+ std::string tmp_name = var_data.varName + "0";
+ tag_name << tmp_name.c_str();
+ }
+ else
+ tag_name << var_data.varName << tstep_num;
+ ErrorCode rval = MB_SUCCESS;
+ tagh = 0;
+ switch (var_data.varDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ rval = mbImpl->tag_get_handle(tag_name.str().c_str(), num_lev, MB_TYPE_OPAQUE, tagh, MB_TAG_DENSE | MB_TAG_CREAT);
+ break;
+ case NC_DOUBLE:
+ case NC_FLOAT:
+ rval = mbImpl->tag_get_handle(tag_name.str().c_str(), num_lev, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT);
+ break;
+ case NC_INT:
+ case NC_SHORT:
+ rval = mbImpl->tag_get_handle(tag_name.str().c_str(), num_lev, MB_TYPE_INTEGER, tagh, MB_TAG_DENSE | MB_TAG_CREAT);
+ break;
+ default:
+ std::cerr << "Unrecognized data type for tag " << tag_name.str() << std::endl;
+ rval = MB_FAILURE;
+ }
+
+ if (MB_SUCCESS == rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.str().c_str());
+
+ return rval;
+}
+
+ErrorCode NCHelper::create_attrib_string(const std::map<std::string, ReadNC::AttData>& attMap, std::string& attVal, std::vector<int>& attLen)
+{
+ int success;
+ std::stringstream ssAtt;
+ unsigned int sz = 0;
+ std::map<std::string, ReadNC::AttData>::const_iterator attIt = attMap.begin();
+ for (; attIt != attMap.end(); ++attIt) {
+ ssAtt << attIt->second.attName;
+ ssAtt << '\0';
+ void* attData = NULL;
+ switch (attIt->second.attDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ sz = attIt->second.attLen;
+ attData = (char *) malloc(sz);
+ success = NCFUNC(get_att_text)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (char*) attData);
+ ERRORS(success, "Failed to read attribute char data.");
+ ssAtt << "char;";
+ break;
+ case NC_DOUBLE:
+ sz = attIt->second.attLen * sizeof(double);
+ attData = (double *) malloc(sz);
+ success = NCFUNC(get_att_double)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (double*) attData);
+ ERRORS(success, "Failed to read attribute double data.");
+ ssAtt << "double;";
+ break;
+ case NC_FLOAT:
+ sz = attIt->second.attLen * sizeof(float);
+ attData = (float *) malloc(sz);
+ success = NCFUNC(get_att_float)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (float*) attData);
+ ERRORS(success, "Failed to read attribute float data.");
+ ssAtt << "float;";
+ break;
+ case NC_INT:
+ sz = attIt->second.attLen * sizeof(int);
+ attData = (int *) malloc(sz);
+ success = NCFUNC(get_att_int)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (int*) attData);
+ ERRORS(success, "Failed to read attribute int data.");
+ ssAtt << "int;";
+ break;
+ case NC_SHORT:
+ sz = attIt->second.attLen * sizeof(short);
+ attData = (short *) malloc(sz);
+ success = NCFUNC(get_att_short)(_fileId, attIt->second.attVarId, attIt->second.attName.c_str(), (short*) attData);
+ ERRORS(success, "Failed to read attribute short data.");
+ ssAtt << "short;";
+ break;
+ default:
+ success = 1;
+ }
+ char* tmpc = (char *) attData;
+ for (unsigned int counter = 0; counter != sz; ++counter)
+ ssAtt << tmpc[counter];
+ free(attData);
+ ssAtt << ';';
+ attLen.push_back(ssAtt.str().size() - 1);
+ }
+ attVal = ssAtt.str();
+
+ return MB_SUCCESS;
+}
+
+void NCHelper::init_dims_with_no_cvars_info()
+{
+ std::vector<std::string>& dimNames = _readNC->dimNames;
+ std::set<std::string>& dummyVarNames = _readNC->dummyVarNames;
+ std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+
+ // Hack: look at all dimensions, and see if we have one that does not appear in the list of varInfo names
+ // right now, candidates are ncol and nbnd
+ // for them, create dummy tags
+ for (unsigned int i = 0; i < dimNames.size(); i++)
+ {
+ // If there is a var with this name, skip, we are fine; if not, create a varInfo...
+ if (varInfo.find(dimNames[i]) != varInfo.end())
+ continue; // We already have a variable with this dimension name
+
+ int sizeTotalVar = varInfo.size();
+ std::string var_name(dimNames[i]);
+ ReadNC::VarData &data = varInfo[var_name];
+ data.varName = std::string(var_name);
+ data.varId =sizeTotalVar;
+ data.varTags.resize(1, 0);
+ data.varDataType = NC_DOUBLE; // Could be int, actually, but we do not really need the type
+ data.varDims.resize(1);
+ data.varDims[0]= (int)i;
+ data.numAtts=0;
+ data.entLoc = ReadNC::ENTLOCSET;
+ dbgOut.tprintf(2, "Dummy varInfo created for dimension %s\n", dimNames[i].c_str());
+ dummyVarNames.insert(dimNames[i]);
+ }
+}
+
+ErrorCode NCHelper::read_variable_to_set_allocate(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums)
+{
+ std::vector<int>& dimVals = _readNC->dimVals;
+ DebugOutput& dbgOut = _readNC->dbgOut;
+
+ ErrorCode rval = MB_SUCCESS;
+
+ for (unsigned int i = 0; i < vdatas.size(); i++) {
+ if ((std::find(vdatas[i].varDims.begin(), vdatas[i].varDims.end(), tDim) != vdatas[i].varDims.end()))
+ vdatas[i].has_t = true;
+
+ for (unsigned int t = 0; t < tstep_nums.size(); t++) {
+ dbgOut.tprintf(2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
+
+ // Get the tag to read into
+ if (!vdatas[i].varTags[t]) {
+ rval = get_tag_to_set(vdatas[i], tstep_nums[t], vdatas[i].varTags[t]);
+ ERRORR(rval, "Trouble getting tag.");
+ }
+
+ // Assume point-based values for now?
+ if (-1 == tDim || dimVals[tDim] <= (int) t)
+ ERRORR(MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number.");
+
+ // Set up the dimensions and counts
+ // First variable dimension is time, if it exists
+ if (vdatas[i].has_t)
+ {
+ if (vdatas[i].varDims.size() != 1)
+ {
+ vdatas[i].readStarts[t].push_back(tstep_nums[t]);
+ vdatas[i].readCounts[t].push_back(1);
+ }
+ else
+ {
+ vdatas[i].readStarts[t].push_back(0);
+ vdatas[i].readCounts[t].push_back(tstep_nums.size());
+ }
+ }
+
+ // Set up other dimensions and counts
+ if (vdatas[i].varDims.empty()) {
+ // Scalar variable
+ vdatas[i].readStarts[t].push_back(0);
+ vdatas[i].readCounts[t].push_back(1);
+ }
+ else {
+ for (unsigned int idx = 0; idx != vdatas[i].varDims.size(); idx++){
+ if (tDim != vdatas[i].varDims[idx]){
+ // Push other variable dimensions, except time, which was already pushed
+ vdatas[i].readStarts[t].push_back(0);
+ vdatas[i].readCounts[t].push_back(dimVals[vdatas[i].varDims[idx]]);
+ }
+ }
+ }
+ std::size_t sz = 1;
+ for (std::size_t idx = 0; idx != vdatas[i].readCounts[t].size(); idx++)
+ sz *= vdatas[i].readCounts[t][idx];
+ vdatas[i].sz = sz;
+ switch (vdatas[i].varDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ vdatas[i].varDatas[t] = new char[sz];
+ break;
+ case NC_DOUBLE:
+ case NC_FLOAT:
+ vdatas[i].varDatas[t] = new double[sz];
+ break;
+ case NC_INT:
+ case NC_SHORT:
+ vdatas[i].varDatas[t] = new int[sz];
+ break;
+ default:
+ std::cerr << "Unrecognized data type for tag " << std::endl;
+ rval = MB_FAILURE;
+ }
+ if (vdatas[i].varDims.size() <= 1)
+ break;
+ }
+ }
+
+ return rval;
+}
+
ErrorCode ScdNCHelper::check_existing_mesh(EntityHandle file_set) {
Interface*& mbImpl = _readNC->mbImpl;
- int (&lDims)[6] = _readNC->lDims;
// Get the number of vertices
int num_verts;
@@ -358,15 +817,8 @@ ErrorCode ScdNCHelper::check_existing_mesh(EntityHandle file_set) {
ErrorCode ScdNCHelper::create_mesh(ScdInterface* scdi, const FileOptions& opts, EntityHandle file_set, Range& faces)
{
Interface*& mbImpl = _readNC->mbImpl;
- 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;
DebugOutput& dbgOut = _readNC->dbgOut;
- int (&locallyPeriodic)[3] = _readNC->locallyPeriodic;
- int (&globallyPeriodic)[3] = _readNC->globallyPeriodic;
ScdParData& parData = _readNC->parData;
Range tmp_range;
@@ -406,7 +858,7 @@ ErrorCode ScdNCHelper::create_mesh(ScdInterface* scdi, const FileOptions& opts,
int di = gDims[3] - gDims[0] + 1;
int dj = gDims[4] - gDims[1] + 1;
assert(dil == (int)ilVals.size() && djl == (int)jlVals.size() &&
- (-1 == lDims[2] || lDims[5]-lDims[2]+1 == (int)klVals.size()));
+ (-1 == lDims[2] || lDims[5]-lDims[2] + 1 == (int)levVals.size()));
#define INDEX(i, j, k) ()
for (kl = lDims[2]; kl <= lDims[5]; kl++) {
k = kl - lDims[2];
@@ -417,7 +869,7 @@ ErrorCode ScdNCHelper::create_mesh(ScdInterface* scdi, const FileOptions& opts,
unsigned int pos = i + j * dil + k * dil * djl;
xc[pos] = ilVals[i];
yc[pos] = jlVals[j];
- zc[pos] = (-1 == lDims[2] ? 0.0 : klVals[k]);
+ zc[pos] = (-1 == lDims[2] ? 0.0 : levVals[k]);
itmp = (!locallyPeriodic[0] && globallyPeriodic[0] && il == gDims[3] ? gDims[0] : il);
*gid_data = (-1 != kl ? kl * di * dj : 0) + jl * di + itmp + 1;
gid_data++;
@@ -468,7 +920,7 @@ ErrorCode ScdNCHelper::read_variables(EntityHandle file_set, std::vector<std::st
ERRORR(rval, "Trouble setting up read variable.");
// Create COORDS tag for quads
- rval = _readNC->create_quad_coordinate_tag(file_set);
+ rval = create_quad_coordinate_tag(file_set);
ERRORR(rval, "Trouble creating coordinate tags to entities quads");
if (!vsetdatas.empty()) {
@@ -488,14 +940,6 @@ ErrorCode ScdNCHelper::read_scd_variable_setup(std::vector<std::string>& var_nam
std::vector<ReadNC::VarData>& vdatas, std::vector<ReadNC::VarData>& vsetdatas)
{
std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;
- int& tMin = _readNC->tMin;
- int& tMax = _readNC->tMax;
- int& iDim = _readNC->iDim;
- int& jDim = _readNC->jDim;
- int& tDim = _readNC->tDim;
- int& iCDim = _readNC->iCDim;
- int& jCDim = _readNC->jCDim;
-
std::map<std::string, ReadNC::VarData>::iterator mit;
// If empty read them all
@@ -536,9 +980,9 @@ ErrorCode ScdNCHelper::read_scd_variable_setup(std::vector<std::string>& var_nam
}
}
- if (tstep_nums.empty() && -1 != tMin) {
+ if (tstep_nums.empty() && nTimeSteps > 0) {
// No timesteps input, get them all
- for (int i = tMin; i <= tMax; i++)
+ for (int i = 0; i < nTimeSteps; i++)
tstep_nums.push_back(i);
}
if (!tstep_nums.empty()) {
@@ -571,16 +1015,9 @@ ErrorCode ScdNCHelper::read_scd_variable_setup(std::vector<std::string>& var_nam
ErrorCode ScdNCHelper::read_scd_variable_to_nonset_allocate(EntityHandle file_set, std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums)
{
Interface*& mbImpl = _readNC->mbImpl;
- std::vector<std::string>& dimNames = _readNC->dimNames;
std::vector<int>& dimVals = _readNC->dimVals;
- int (&lDims)[6] = _readNC->lDims;
- int (&lCDims)[6] = _readNC->lCDims;
- int& tDim = _readNC->tDim;
DebugOutput& dbgOut = _readNC->dbgOut;
bool& isParallel = _readNC->isParallel;
- #ifdef USE_MPI
- ParallelComm*& myPcomm = _readNC->myPcomm;
-#endif
ErrorCode rval = MB_SUCCESS;
@@ -606,8 +1043,8 @@ ErrorCode ScdNCHelper::read_scd_variable_to_nonset_allocate(EntityHandle file_se
#ifdef USE_MPI
moab::Range faces_owned;
- if (isParallel)
- {
+ if (isParallel) {
+ ParallelComm*& myPcomm = _readNC->myPcomm;
rval = myPcomm->filter_pstatus(faces, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &faces_owned);
ERRORR(rval, "Trouble getting owned faces in set.");
}
@@ -616,24 +1053,14 @@ ErrorCode ScdNCHelper::read_scd_variable_to_nonset_allocate(EntityHandle file_se
#endif
for (unsigned int i = 0; i < vdatas.size(); i++) {
+ vdatas[i].numLev = nLevels;
+
for (unsigned int t = 0; t < tstep_nums.size(); t++) {
dbgOut.tprintf(2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
- std::vector<std::string>::iterator vit;
- int idx_lev = -1;
- int idx_ilev = -1;
- if ((vit = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end())
- idx_lev = vit - dimNames.begin();
- if ((vit = std::find(dimNames.begin(), dimNames.end(), "ilev")) != dimNames.end())
- idx_ilev = vit - dimNames.begin();
- if (std::find(vdatas[i].varDims.begin(), vdatas[i].varDims.end(), idx_lev) != vdatas[i].varDims.end())
- vdatas[i].numLev = dimVals[idx_lev];
- else if (std::find(vdatas[i].varDims.begin(), vdatas[i].varDims.end(), idx_ilev) != vdatas[i].varDims.end())
- vdatas[i].numLev = dimVals[idx_ilev];
-
// Get the tag to read into
if (!vdatas[i].varTags[t]) {
- rval = _readNC->get_tag_to_nonset(vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev);
+ rval = get_tag_to_nonset(vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev);
ERRORR(rval, "Trouble getting tag.");
}
@@ -838,6 +1265,66 @@ ErrorCode ScdNCHelper::read_scd_variable_to_nonset(EntityHandle file_set, std::v
return rval;
}
+ErrorCode ScdNCHelper::create_quad_coordinate_tag(EntityHandle file_set) {
+ Interface*& mbImpl = _readNC->mbImpl;
+ bool& isParallel = _readNC->isParallel;
+
+ Range ents;
+ ErrorCode rval = mbImpl->get_entities_by_type(file_set, moab::MBQUAD, ents);
+ ERRORR(rval, "Trouble getting QUAD entity.");
+
+ std::size_t numOwnedEnts = 0;
+#ifdef USE_MPI
+ Range ents_owned;
+ if (isParallel) {
+ ParallelComm*& myPcomm = _readNC->myPcomm;
+ rval = myPcomm->filter_pstatus(ents, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &ents_owned);
+ ERRORR(rval, "Trouble getting owned QUAD entity.");
+ numOwnedEnts = ents_owned.size();
+ }
+ else
+ {
+ numOwnedEnts = ents.size();
+ ents_owned = ents;
+ }
+#else
+ numOwnedEnts = ents.size();
+#endif
+
+ if (numOwnedEnts == 0)
+ return MB_SUCCESS;
+
+ assert(numOwnedEnts == ilCVals.size() * jlCVals.size());
+ std::vector<double> coords(numOwnedEnts * 3);
+ std::size_t pos = 0;
+ for (std::size_t j = 0; j != jlCVals.size(); ++j) {
+ for (std::size_t i = 0; i != ilCVals.size(); ++i) {
+ pos = j * ilCVals.size() * 3 + i * 3;
+ coords[pos] = ilCVals[i];
+ coords[pos + 1] = jlCVals[j];
+ coords[pos + 2] = 0.0;
+ }
+ }
+ std::string tag_name = "COORDS";
+ Tag tagh = 0;
+ rval = mbImpl->tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating COORDS tag.");
+
+ void *data;
+ int count;
+#ifdef USE_MPI
+ rval = mbImpl->tag_iterate(tagh, ents_owned.begin(), ents_owned.end(), count, data);
+#else
+ rval = mbImpl->tag_iterate(tagh, ents.begin(), ents.end(), count, data);
+#endif
+ ERRORR(rval, "Failed to get COORDS tag iterator.");
+ assert(count == (int)numOwnedEnts);
+ double* quad_data = (double*) data;
+ std::copy(coords.begin(), coords.end(), quad_data);
+
+ return MB_SUCCESS;
+}
+
ErrorCode UcdNCHelper::read_variables(EntityHandle file_set, std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
{
std::vector<ReadNC::VarData> vdatas;
diff --git a/src/io/NCHelper.hpp b/src/io/NCHelper.hpp
index f59663e..83fcec1 100644
--- a/src/io/NCHelper.hpp
+++ b/src/io/NCHelper.hpp
@@ -17,7 +17,8 @@ namespace moab {
class NCHelper
{
public:
- NCHelper(ReadNC* readNC, int fileId) : _readNC(readNC), _fileId(fileId) {}
+ NCHelper(ReadNC* readNC, int fileId) :_readNC(readNC), _fileId(fileId),
+ nTimeSteps(0), nLevels(1), tDim(-1), levDim(-1) {}
virtual ~NCHelper() {}
//! Get appropriate helper instance for ReadNC class
@@ -30,6 +31,10 @@ public:
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;
+ //! Create NC conventional tags
+ ErrorCode create_conventional_tags(ScdInterface* scdi, EntityHandle file_set,
+ const std::vector<int>& tstep_nums);
+
protected:
//! Read set variables, common to scd mesh and ucd mesh
ErrorCode read_variable_to_set(EntityHandle file_set, std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums);
@@ -37,20 +42,62 @@ protected:
//! Convert variables in place
ErrorCode convert_variable(ReadNC::VarData& var_data, int tstep_num);
+ ErrorCode read_coordinate(const char* var_name, int lmin, int lmax,
+ std::vector<double>& cvals);
+
+ ErrorCode get_tag_to_set(ReadNC::VarData& var_data, int tstep_num, Tag& tagh);
+
+ ErrorCode get_tag_to_nonset(ReadNC::VarData& var_data, int tstep_num, Tag& tagh, int num_lev);
+
+ //! 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
+ //! the next'. attLen stores the end position for each name/data
+ //! type/ value.
+ ErrorCode create_attrib_string(const std::map<std::string, ReadNC::AttData>& attMap,
+ std::string& attString,
+ std::vector<int>& attLen);
+
+ //! Init info for dimensions that don't have corresponding
+ //! coordinate variables - this info is used for creating tags
+ void init_dims_with_no_cvars_info();
+
private:
//! Used by read_variable_to_set()
ErrorCode read_variable_to_set_allocate(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums);
protected:
ReadNC* _readNC;
+
int _fileId;
+
+ //! Dimensions of time and level
+ int nTimeSteps, nLevels;
+
+ //! Values for time and level
+ std::vector<double> tVals, levVals;
+
+ //! Dimension numbers for time and level
+ int tDim, levDim;
};
//! Child helper class for scd mesh, e.g. CAM_EL or CAM_FV
class ScdNCHelper : public NCHelper
{
public:
- ScdNCHelper(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+ ScdNCHelper(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId),
+ iDim(-1), jDim(-1), iCDim(-1), jCDim(-1)
+ {
+ for (unsigned int i = 0; i < 6; i++) {
+ gDims[i] = -1;
+ lDims[i] = -1;
+ gCDims[i] = -1;
+ lCDims[i] = -1;
+ }
+
+ locallyPeriodic[0] = locallyPeriodic[1] = 0;
+ globallyPeriodic[0] = globallyPeriodic[1] = 0;
+ }
virtual ~ScdNCHelper() {}
private:
@@ -73,6 +120,9 @@ private:
ErrorCode read_scd_variable_to_nonset(EntityHandle file_set, std::vector<ReadNC::VarData>& vdatas,
std::vector<int>& tstep_nums);
+ //! Create COORDS tag for quads coordinate
+ ErrorCode create_quad_coordinate_tag(EntityHandle file_set);
+
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;
@@ -83,6 +133,37 @@ private:
tmp_data[j*nik + i*nk + k] = source[k*nij + j*ni + i];
return MB_SUCCESS;
}
+
+protected:
+ //! Dimensions of global grid in file
+ int gDims[6];
+
+ //! Dimensions of my local part of grid
+ int lDims[6];
+
+ //! Center dimensions of global grid in file
+ int gCDims[6];
+
+ //! Center dimensions of my local part of grid
+ int lCDims[6];
+
+ //! Values for i/j
+ std::vector<double> ilVals, jlVals;
+
+ //! Center values for i/j
+ std::vector<double> ilCVals, jlCVals;
+
+ //! Dimension numbers for i/j
+ int iDim, jDim;
+
+ //! Center dimension numbers for i/j
+ int iCDim, jCDim;
+
+ //! Whether mesh is locally periodic in i or j
+ int locallyPeriodic[2];
+
+ //! Whether mesh is globally periodic in i or j
+ int globallyPeriodic[2];
};
//! Child helper class for ucd mesh, e.g. CAM_SE (HOMME) or MPAS
@@ -90,9 +171,9 @@ class UcdNCHelper : public NCHelper
{
public:
UcdNCHelper(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId),
- cDim(-1), eDim(-1), vDim(-1), levDim(-1),
nCells(0), nEdges(0), nVertices(0),
- nLocalCells(0), nLocalEdges(0), nLocalVertices(0) {}
+ nLocalCells(0), nLocalEdges(0), nLocalVertices(0),
+ cDim(-1), eDim(-1), vDim(-1) {}
virtual ~UcdNCHelper() {}
private:
@@ -139,20 +220,23 @@ protected:
return MB_SUCCESS;
}
- //! Dimension numbers for nCells, nEdges, nVertices, nLevels
- int cDim, eDim, vDim, levDim;
-
- //! Coordinate values for vertices
- std::vector<double> xVertVals, yVertVals, zVertVals;
-
+ //! Dimensions of global grid in file
int nCells;
int nEdges;
int nVertices;
+ //! Dimensions of my local part of grid
int nLocalCells;
int nLocalEdges;
int nLocalVertices;
+ //! Coordinate values for vertices
+ std::vector<double> xVertVals, yVertVals, zVertVals;
+
+ //! Dimension numbers for nCells, nEdges and nVertices
+ int cDim, eDim, vDim;
+
+ //! Local global ID for cells, edges and vertices
Range localGidCells, localGidEdges, localGidVerts;
};
diff --git a/src/io/NCHelperEuler.cpp b/src/io/NCHelperEuler.cpp
index 8de0818..51094dd 100644
--- a/src/io/NCHelperEuler.cpp
+++ b/src/io/NCHelperEuler.cpp
@@ -54,78 +54,51 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
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)[3] = _readNC->locallyPeriodic;
- int (&globallyPeriodic)[3] = _readNC->globallyPeriodic;
ScdParData& parData = _readNC->parData;
-#ifdef USE_MPI
- ParallelComm*& myPcomm = _readNC->myPcomm;
-#endif
- // look for names of center i/j dimensions
+ // Look for names of center i/j dimensions
+ // First i
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())
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lon")) != dimNames.end())
idx = vit - dimNames.begin();
else {
- ERRORR(MB_FAILURE, "Couldn't find center i variable.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lon' dimension.");
}
iCDim = idx;
+ gCDims[0] = 0;
+ gCDims[3] = dimVals[idx] - 1;
- // 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)
+ // Check i periodicity and set globallyPeriodic[0]
+ std::vector<double> til_vals(2);
+ ErrorCode rval = read_coordinate("lon", gCDims[3] - 1, gCDims[3], til_vals);
+ ERRORR(rval, "Trouble reading 'lon' variable.");
+ if (std::fabs(2 * til_vals[1] - til_vals[0] - 360) < 0.001)
globallyPeriodic[0] = 1;
- // now we can set gCDims and gDims for i
- gCDims[0] = 0;
+ // Now we can set gDims for i
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
+ 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())
+ // Then j
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lat")) != dimNames.end())
idx = vit - dimNames.begin();
else {
- ERRORR(MB_FAILURE, "Couldn't find center j variable.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lat' dimension.");
}
jCDim = idx;
-
- // for Eul models, will always be non-periodic in j
gCDims[1] = 0;
- gDims[1] = 0;
gCDims[4] = dimVals[idx] - 1;
+
+ // For Eul models, will always be non-periodic in j
+ gDims[1] = 0;
gDims[4] = gCDims[4] + 1;
- // try a truly 2d mesh
+ // Try a truly 2D mesh
gDims[2] = -1;
gDims[5] = -1;
@@ -135,16 +108,32 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
idx = vit - dimNames.begin();
else {
- ERRORR(MB_FAILURE, "Couldn't find time dimension.");
+ ERRORR(MB_FAILURE, "Couldn't find 'time' or 't' dimension.");
}
tDim = idx;
- tMax = dimVals[idx] - 1;
- tMin = 0;
- tName = dimNames[idx];
+ nTimeSteps = dimVals[idx];
- // parse options to get subset
- if (isParallel) {
+ // Look for level dimension
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else if ((vit = std::find(dimNames.begin(), dimNames.end(), "ilev")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension.");
+ }
+ levDim = idx;
+ nLevels = dimVals[idx];
+
+ // Parse options to get subset
+ int rank = 0, procs = 1;
#ifdef USE_MPI
+ if (isParallel) {
+ ParallelComm*& myPcomm = _readNC->myPcomm;
+ rank = myPcomm->proc_config().proc_rank();
+ procs = myPcomm->proc_config().proc_size();
+ }
+#endif
+ if (procs > 1) {
for (int i = 0; i < 6; i++)
parData.gDims[i] = gDims[i];
for (int i = 0; i < 3; i++)
@@ -152,9 +141,7 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
parData.partMethod = partMethod;
int pdims[3];
- rval = ScdInterface::compute_partition(myPcomm->proc_config().proc_size(),
- myPcomm->proc_config().proc_rank(),
- parData, lDims, locallyPeriodic, pdims);
+ rval = ScdInterface::compute_partition(procs, rank, parData, lDims, locallyPeriodic, pdims);
if (MB_SUCCESS != rval)
return rval;
for (int i = 0; i < 3; i++)
@@ -163,9 +150,8 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
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)
+ if (0 == rank)
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++)
@@ -178,92 +164,83 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
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);
+ // Now get actual coordinate values for vertices and cell centers
+ lCDims[0] = lDims[0];
+ if (locallyPeriodic[0])
+ // If locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem coords
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;
- }
- }
+ else
+ lCDims[3] = lDims[3] - 1;
- lCDims[0] = lDims[0];
- lCDims[4] = lDims[4] - 1;
+ // For Eul models, will always be non-periodic in j
lCDims[1] = lDims[1];
+ lCDims[4] = lDims[4] - 1;
- if (-1 != lDims[1]) {
+ // Resize vectors to store values later
+ if (-1 != lDims[0])
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ if (-1 != lCDims[0])
+ ilCVals.resize(lCDims[3] - lCDims[0] + 1);
+ if (-1 != lDims[1])
jlVals.resize(lDims[4] - lDims[1] + 1);
+ if (-1 != lCDims[1])
jlCVals.resize(lCDims[4] - lCDims[1] + 1);
- }
+ if (nTimeSteps > 0)
+ tVals.resize(nTimeSteps);
- if (-1 != tMin)
- tVals.resize(tMax - tMin + 1);
-
- // now read coord values
+ // 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.");
+ if (-1 != lCDims[0]) {
+ if ((vmit = varInfo.find("lon")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = read_coordinate("lon", lCDims[0], lCDims[3], ilCVals);
+ ERRORR(rval, "Trouble reading 'lon' variable.");
}
else {
- ERRORR(MB_FAILURE, "Couldn't find lon coordinate.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lon' variable.");
}
}
- 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.");
+ if (-1 != lCDims[1]) {
+ if ((vmit = varInfo.find("lat")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = read_coordinate("lat", lCDims[1], lCDims[4], jlCVals);
+ ERRORR(rval, "Trouble reading 'lat' variable.");
}
else {
- ERRORR(MB_FAILURE, "Couldn't find lat coordinate.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lat' variable.");
}
}
- if (lDims[0] != -1) {
- if ((vmit = varInfo.find(iCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ if (-1 != lDims[0]) {
+ if ((vmit = varInfo.find("lon")) != 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
+ // 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.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lon' variable.");
}
}
- if (lDims[1] != -1) {
- if ((vmit = varInfo.find(jCName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ if (-1 != lDims[1]) {
+ if ((vmit = varInfo.find("lat")) != 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
+ rval = 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;
+ std::size_t i = 0;
double gwSum = -1;
for (i = 1; i != gwVals.size() + 1; i++) {
- gwSum = gwSum + gwVals[i - 1];
+ gwSum += gwVals[i - 1];
jlVals[i] = std::asin(gwSum) * 180 / M_PI;
}
- jlVals[i] = 90.0; // using value of i after loop exits.
+ jlVals[i] = 90.0; // Using value of i after loop exits.
}
else {
std::string gwName("gw");
@@ -272,66 +249,68 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
// 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
+ rval = 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];
+ gwSum += gwVals[i - 1];
jlVals[i] = std::asin(gwSum) * 180 / M_PI;
}
}
- // or if it's the last row
+ // 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.");
+ rval = 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];
- }
+ for (int j = 0; j != lDims[1] - 1; j++)
+ gwSum += gwVals[j];
std::size_t i = 0;
for (; i != jlVals.size() - 1; i++) {
- gwSum = gwSum + gwVals[lDims[1] - 1 + i];
+ 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.
+ jlVals[i] = 90.0; // Using value of i after loop exits.
}
- // it's in the middle
+ // 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.");
+ rval = 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];
- }
+ for (int j = 0; j != start - 1; j++)
+ gwSum += gwVals[j];
std::size_t i = 0;
for (; i != jlVals.size(); i++) {
- gwSum = gwSum + gwVals[start - 1 + i];
+ gwSum += gwVals[start - 1 + i];
jlVals[i] = std::asin(gwSum) * 180 / M_PI;
}
}
}
}
else {
- ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lat' variable.");
}
}
// Store time coordinate values in tVals
- 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.");
+ if (nTimeSteps > 0) {
+ if ((vmit = varInfo.find("time")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = read_coordinate("time", 0, nTimeSteps - 1, tVals);
+ ERRORR(rval, "Trouble reading 'time' variable.");
+ }
+ else if ((vmit = varInfo.find("t")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = read_coordinate("t", 0, nTimeSteps - 1, tVals);
+ ERRORR(rval, "Trouble reading 't' variable.");
}
else {
// If expected time variable is not available, set dummy time coordinate values to tVals
- for (int t = tMin; t <= tMax; t++)
+ for (int t = 0; t < nTimeSteps; t++)
tVals.push_back((double)t);
}
}
@@ -340,7 +319,7 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
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
+ // 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;
@@ -360,7 +339,7 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
int val_len = 0;
for (unsigned int i = 0; i != ijdimNames.size(); i++) {
tag_name = ijdimNames[i];
- void * val = NULL;
+ void* val = NULL;
if (tag_name == "__slon") {
val = &ilVals[0];
val_len = ilVals.size();
@@ -380,7 +359,7 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
Tag tagh = 0;
DataType data_type;
- // assume all has same data type as lon
+ // Assume all has same data type as lon
switch (varInfo["lon"].varDataType) {
case NC_BYTE:
case NC_CHAR:
@@ -469,9 +448,8 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
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();
+ // Hack: create dummy tags, if needed, for variables with no corresponding variables
+ init_dims_with_no_cvars_info();
return MB_SUCCESS;
}
diff --git a/src/io/NCHelperFV.cpp b/src/io/NCHelperFV.cpp
index 8b9be3f..ccae4e3 100644
--- a/src/io/NCHelperFV.cpp
+++ b/src/io/NCHelperFV.cpp
@@ -47,76 +47,50 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
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)[3] = _readNC->locallyPeriodic;
- int (&globallyPeriodic)[3] = _readNC->globallyPeriodic;
ScdParData& parData = _readNC->parData;
-#ifdef USE_MPI
- ParallelComm*& myPcomm = _readNC->myPcomm;
-#endif
+ // Look for names of i/j dimensions
+ // First i
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.");
+ ERRORR(MB_FAILURE, "Couldn't find 'slon' variable.");
}
-
iDim = idx;
- gDims[3] = dimVals[idx] - 1;
gDims[0] = 0;
- iName = dimNames[idx];
+ gDims[3] = dimVals[idx] - 1;
+ // Then j
if ((vit = std::find(dimNames.begin(), dimNames.end(), "slat")) != dimNames.end())
idx = vit - dimNames.begin();
else {
- ERRORR(MB_FAILURE, "Couldn't find slat variable.");
+ 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];
+ gDims[4] = dimVals[idx] - 1 + 2; // Add 2 for the pole points
- // look for names of center i/j dimensions
+ // Look for names of center i/j dimensions
+ // First i
if ((vit = std::find(dimNames.begin(), dimNames.end(), "lon")) != dimNames.end())
idx = vit - dimNames.begin();
else {
- ERRORR(MB_FAILURE, "Couldn't find lon variable.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lon' variable.");
}
iCDim = idx;
- gCDims[3] = dimVals[idx] - 1;
gCDims[0] = 0;
- iCName = dimNames[idx];
+ gCDims[3] = dimVals[idx] - 1;
- // check and set globallyPeriodic[0]
+ // Check i periodicity 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.");
+ ErrorCode rval = read_coordinate("lon", gCDims[3] - 1, gCDims[3], til_vals);
+ ERRORR(rval, "Trouble reading 'lon' variable.");
if (std::fabs(2 * til_vals[1] - til_vals[0] - 360) < 0.001)
globallyPeriodic[0] = 1;
if (globallyPeriodic[0])
@@ -124,39 +98,55 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
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
-
+ // Then j
if ((vit = std::find(dimNames.begin(), dimNames.end(), "lat")) != dimNames.end())
idx = vit - dimNames.begin();
else {
- ERRORR(MB_FAILURE, "Couldn't find lat variable.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lat' dimension.");
}
jCDim = idx;
- gCDims[4] = dimVals[idx] - 1;
gCDims[1] = 0;
- jCName = dimNames[idx];
+ gCDims[4] = dimVals[idx] - 1;
+
+ // For FV models, will always be non-periodic in j
+ assert(gDims[4] == gCDims[4] + 1);
+
+ // Try a truly 2D mesh
+ gDims[2] = -1;
+ gDims[5] = -1;
// Look for time dimension
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 dimension.");
+ ERRORR(MB_FAILURE, "Couldn't find 'time' or 't' dimension.");
}
tDim = idx;
- tMax = dimVals[idx] - 1;
- tMin = 0;
- tName = dimNames[idx];
+ nTimeSteps = dimVals[idx];
- // parse options to get subset
- if (isParallel) {
+ // Get number of levels
+ if ((vit = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else if ((vit = std::find(dimNames.begin(), dimNames.end(), "ilev")) != dimNames.end())
+ idx = vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension.");
+ }
+ levDim = idx;
+ nLevels = dimVals[idx];
+
+ // Parse options to get subset
+ int rank = 0, procs = 1;
#ifdef USE_MPI
+ if (isParallel) {
+ ParallelComm*& myPcomm = _readNC->myPcomm;
+ rank = myPcomm->proc_config().proc_rank();
+ procs = myPcomm->proc_config().proc_size();
+ }
+#endif
+ if (procs > 1) {
for (int i = 0; i < 6; i++)
parData.gDims[i] = gDims[i];
for (int i = 0; i < 3; i++)
@@ -164,9 +154,7 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
parData.partMethod = partMethod;
int pdims[3];
- rval = ScdInterface::compute_partition(myPcomm->proc_config().proc_size(),
- myPcomm->proc_config().proc_rank(),
- parData, lDims, locallyPeriodic, pdims);
+ rval = ScdInterface::compute_partition(procs, rank, parData, lDims, locallyPeriodic, pdims);
if (MB_SUCCESS != rval)
return rval;
for (int i = 0; i < 3; i++)
@@ -175,154 +163,150 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
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)
+ if (0 == rank)
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);
+ // Now get actual coordinate values for vertices and cell centers
+ lCDims[0] = lDims[0];
+ if (locallyPeriodic[0])
+ // If locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem coords
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;
- }
- }
+ else
+ lCDims[3] = lDims[3] - 1;
- lCDims[0] = lDims[0];
- lCDims[4] = lDims[4] - 1;
+ // For FV models, will always be non-periodic in j
lCDims[1] = lDims[1];
+ lCDims[4] = lDims[4] - 1;
- if (-1 != lDims[1]) {
+ // Resize vectors to store values later
+ if (-1 != lDims[0])
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ if (-1 != lCDims[0])
+ ilCVals.resize(lCDims[3] - lCDims[0] + 1);
+ if (-1 != lDims[1])
jlVals.resize(lDims[4] - lDims[1] + 1);
+ if (-1 != lCDims[1])
jlCVals.resize(lCDims[4] - lCDims[1] + 1);
- }
-
- if (-1 != tMin)
- tVals.resize(tMax - tMin + 1);
+ if (nTimeSteps > 0)
+ tVals.resize(nTimeSteps);
- // ... then read actual values
+ // Now read coord 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.");
+ if (-1 != lCDims[0]) {
+ if ((vmit = varInfo.find("lon")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = read_coordinate("lon", lCDims[0], lCDims[3], ilCVals);
+ ERRORR(rval, "Trouble reading 'lon' variable.");
}
else {
- ERRORR(MB_FAILURE, "Couldn't find lon coordinate.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lon' variable.");
}
}
- 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.");
+ if (-1 != lCDims[1]) {
+ if ((vmit = varInfo.find("lat")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = read_coordinate("lat", lCDims[1], lCDims[4], jlCVals);
+ ERRORR(rval, "Trouble reading 'lat' variable.");
}
else {
- ERRORR(MB_FAILURE, "Couldn't find lat coordinate.");
+ ERRORR(MB_FAILURE, "Couldn't find 'lat' variable.");
}
}
- if (lDims[0] != -1) {
- if ((vmit = varInfo.find(iName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
- // last column
+ if (-1 != lDims[0]) {
+ if ((vmit = varInfo.find("slon")) != 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];
+ assert(lDims[3] == gDims[3] + 1);
+ std::vector<double> dummyVar(lDims[3] - lDims[0]);
+ rval = read_coordinate("slon", lDims[0], lDims[3] - 1, dummyVar);
+ double dif = dummyVar[1] - dummyVar[0];
std::size_t i;
- for (i = 0; i != til_vals.size(); i++)
- ilVals[i] = til_vals[i];
+ for (i = 0; i != dummyVar.size(); i++)
+ ilVals[i] = dummyVar[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.");
+ rval = read_coordinate("slon", lDims[0], lDims[3], ilVals);
+ ERRORR(rval, "Trouble reading 'slon' variable.");
}
}
else {
- ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ ERRORR(MB_FAILURE, "Couldn't find 'slon' variable.");
}
}
- if (lDims[1] != -1) {
- if ((vmit = varInfo.find(jName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ if (-1 != lDims[1]) {
+ if ((vmit = varInfo.find("slat")) != 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
+ rval = read_coordinate("slat", lDims[1], lDims[4] - 2, dummyVar);
+ ERRORR(rval, "Trouble reading 'slat' variable.");
+ // Copy the correct piece
jlVals[0] = -90.0;
- unsigned int i = 0;
+ std::size_t 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.
+ 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
+ // 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
+ rval = read_coordinate("slat", lDims[1], lDims[4] - 1, dummyVar);
+ ERRORR(rval, "Trouble reading 'slat' 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
+ // 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
+ rval = read_coordinate("slat", lDims[1] - 1, lDims[4] - 2, dummyVar);
+ ERRORR(rval, "Trouble reading 'slat' 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.
+ jlVals[i] = 90.0; // Using value of i after loop exits.
}
- // it's in the middle
+ // 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.");
+ rval = read_coordinate("slat", lDims[1] - 1, lDims[4] - 1, jlVals);
+ ERRORR(rval, "Trouble reading 'slat' variable.");
}
}
}
else {
- ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ ERRORR(MB_FAILURE, "Couldn't find 'slat' variable.");
}
}
// Store time coordinate values in tVals
- 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.");
+ if (nTimeSteps > 0) {
+ if ((vmit = varInfo.find("time")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = read_coordinate("time", 0, nTimeSteps - 1, tVals);
+ ERRORR(rval, "Trouble reading 'time' variable.");
+ }
+ else if ((vmit = varInfo.find("t")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
+ rval = read_coordinate("t", 0, nTimeSteps - 1, tVals);
+ ERRORR(rval, "Trouble reading 't' variable.");
}
else {
// If expected time variable is not available, set dummy time coordinate values to tVals
- for (int t = tMin; t <= tMax; t++)
+ for (int t = 0; t < nTimeSteps; t++)
tVals.push_back((double)t);
}
}
@@ -331,7 +315,7 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
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
+ // 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;
@@ -357,7 +341,7 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
int val_len = 0;
for (unsigned int i = 0; i != ijdimNames.size(); i++) {
tag_name = ijdimNames[i];
- void * val = NULL;
+ void* val = NULL;
if (tag_name == "__slon") {
val = &ilVals[0];
val_len = ilVals.size();
@@ -377,7 +361,7 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
Tag tagh = 0;
DataType data_type;
- // assume all has same data type as lon
+ // Assume all has same data type as lon
switch (varInfo["lon"].varDataType) {
case NC_BYTE:
case NC_CHAR:
@@ -466,9 +450,8 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
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();
+ // Hack: create dummy tags, if needed, for dimensions with no corresponding variables
+ init_dims_with_no_cvars_info();
return MB_SUCCESS;
}
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