[MOAB-dev] commit/MOAB: danwu: Merged master into error_handling_enhancement

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Fri Jun 6 14:01:57 CDT 2014


1 new commit in MOAB:

https://bitbucket.org/fathomteam/moab/commits/d7b65921e6db/
Changeset:   d7b65921e6db
Branch:      error_handling_enhancement
User:        danwu
Date:        2014-06-06 21:01:52
Summary:     Merged master into error_handling_enhancement
Affected #:  4 files

diff --git a/src/io/NCHelperGCRM.cpp b/src/io/NCHelperGCRM.cpp
index 25fe9ac..f673aff 100644
--- a/src/io/NCHelperGCRM.cpp
+++ b/src/io/NCHelperGCRM.cpp
@@ -18,12 +18,11 @@
 
 namespace moab {
 
-const int DEFAULT_MAX_EDGES_PER_CELL = 6;
+// GCRM cells are either pentagons or hexagons, and pentagons are always padded to hexagons
+const int EDGES_PER_CELL = 6;
 
 NCHelperGCRM::NCHelperGCRM(ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
 : UcdNCHelper(readNC, fileId, opts, fileSet)
-, maxEdgesPerCell(DEFAULT_MAX_EDGES_PER_CELL)
-, numCellGroups(0)
 , createGatherSet(false)
 {
   // Ignore variables containing topological information
@@ -38,7 +37,7 @@ bool NCHelperGCRM::can_read_file(ReadNC* readNC)
 {
   std::vector<std::string>& dimNames = readNC->dimNames;
 
-  // If dimension name "vertexDegree" exists then it should be the GCRM grid
+  // If dimension name "cells" exists then it should be the GCRM grid
   if (std::find(dimNames.begin(), dimNames.end(), std::string("cells")) != dimNames.end())
     return true;
 
@@ -94,7 +93,7 @@ ErrorCode NCHelperGCRM::init_mesh_vals()
   vDim = idx;
   nVertices = dimLens[idx];
 
-  // Get number of vertex levels
+  // Get number of layers
   if ((vit = std::find(dimNames.begin(), dimNames.end(), "layers")) != dimNames.end())
     idx = vit - dimNames.begin();
   else {
@@ -179,14 +178,6 @@ ErrorCode NCHelperGCRM::check_existing_mesh()
   if (noMesh) {
     ErrorCode rval;
 
-    // Restore numCellGroups
-    if (0 == numCellGroups) {
-      Tag numCellGroupsTag;
-      rval = mbImpl->tag_get_handle("__NUM_CELL_GROUPS", 1, MB_TYPE_INTEGER, numCellGroupsTag);
-      if (MB_SUCCESS == rval)
-        rval = mbImpl->tag_get_data(numCellGroupsTag, &_fileSet, 1, &numCellGroups);
-    }
-
     if (localGidVerts.empty()) {
       // Get all vertices from tmp_set (it is the input set in no_mesh scenario)
       Range local_verts;
@@ -247,14 +238,6 @@ ErrorCode NCHelperGCRM::check_existing_mesh()
         // Restore localGidCells
         std::copy(gids.rbegin(), gids.rend(), range_inserter(localGidCells));
         nLocalCells = localGidCells.size();
-
-        if (numCellGroups > 1) {
-          // Restore cellHandleToGlobalID map
-          Range::const_iterator rit;
-          int i;
-          for (rit = local_cells.begin(), i = 0; rit != local_cells.end(); ++rit, i++)
-            cellHandleToGlobalID[*rit] = gids[i];
-        }
       }
     }
   }
@@ -264,9 +247,7 @@ ErrorCode NCHelperGCRM::check_existing_mesh()
 
 ErrorCode NCHelperGCRM::create_mesh(Range& faces)
 {
-  Interface*& mbImpl = _readNC->mbImpl;
   int& gatherSetRank = _readNC->gatherSetRank;
-  bool& noMixedElements = _readNC->noMixedElements;
   bool& noEdges = _readNC->noEdges;
   DebugOutput& dbgOut = _readNC->dbgOut;
 
@@ -317,9 +298,8 @@ ErrorCode NCHelperGCRM::create_mesh(Range& faces)
   int verticesOnCellVarId;
   int success = NCFUNC(inq_varid)(_fileId, "cell_corners", &verticesOnCellVarId);
   ERRORS(success, "Failed to get variable id of cell_corners.");
-  std::vector<int> vertices_on_local_cells(nLocalCells * maxEdgesPerCell);
+  std::vector<int> vertices_on_local_cells(nLocalCells * EDGES_PER_CELL);
   dbgOut.tprintf(1, " nLocalCells = %d\n", (int)nLocalCells);
-  dbgOut.tprintf(1, " maxEdgesPerCell = %d\n", (int)maxEdgesPerCell);
   dbgOut.tprintf(1, " vertices_on_local_cells.size() = %d\n", (int)vertices_on_local_cells.size());
 #ifdef PNETCDF_FILE
   size_t nb_reads = localGidCells.psize();
@@ -337,7 +317,7 @@ ErrorCode NCHelperGCRM::create_mesh(Range& faces)
     dbgOut.tprintf(1, " cell_corners   endh = %d\n", (int)endh);
     NCDF_SIZE read_starts[2] = {static_cast<NCDF_SIZE>(starth - 1), 0};
     NCDF_SIZE read_counts[2] = {static_cast<NCDF_SIZE>(endh - starth + 1), 
-                                static_cast<NCDF_SIZE>(maxEdgesPerCell)};
+                                static_cast<NCDF_SIZE>(EDGES_PER_CELL)};
 
     // Do a partial read in each subrange
 #ifdef PNETCDF_FILE
@@ -350,7 +330,7 @@ ErrorCode NCHelperGCRM::create_mesh(Range& faces)
     ERRORS(success, "Failed to read cell_corners data in a loop");
 
     // Increment the index for next subrange
-    indexInArray += (endh - starth + 1) * maxEdgesPerCell;
+    indexInArray += (endh - starth + 1) * EDGES_PER_CELL;
   }
 
 #ifdef PNETCDF_FILE
@@ -359,11 +339,25 @@ ErrorCode NCHelperGCRM::create_mesh(Range& faces)
   ERRORS(success, "Failed on wait_all.");
 #endif
 
-  // GCRM is 0 based, convert vertex indices from 0 to 1 based
-  for (std::size_t idx = 0; idx < vertices_on_local_cells.size(); idx++) {
-      vertices_on_local_cells[idx] += 1;
+  // Correct vertices_on_local_cells array. Pentagons as hexagons should have
+  // a connectivity like 123455 and not 122345
+  for (int local_cell_idx = 0; local_cell_idx < nLocalCells; local_cell_idx++) {
+    int* pvertex = &vertices_on_local_cells[local_cell_idx * EDGES_PER_CELL];
+    for (int k = 0; k < EDGES_PER_CELL - 2; k++) {
+      if (*(pvertex + k) == *(pvertex + k + 1)) {
+        // Shift the connectivity
+        for (int kk = k + 1; kk < EDGES_PER_CELL - 1; kk++)
+          *(pvertex + kk) = *(pvertex + kk + 1);
+        // No need to try next k
+        break;
+      }
+    }
   }
 
+  // GCRM is 0 based, convert vertex indices from 0 to 1 based
+  for (std::size_t idx = 0; idx < vertices_on_local_cells.size(); idx++)
+    vertices_on_local_cells[idx] += 1;
+
   // Create local vertices
   EntityHandle start_vertex;
   ErrorCode rval = create_local_vertices(vertices_on_local_cells, start_vertex);
@@ -375,22 +369,9 @@ ErrorCode NCHelperGCRM::create_mesh(Range& faces)
     ERRORR(rval, "Failed to create local edges for GCRM mesh.");
   }
 
-  // Create local cells, either unpadded or padded
-  if (noMixedElements) {
-    rval = create_padded_local_cells(vertices_on_local_cells, start_vertex, faces);
-    ERRORR(rval, "Failed to create padded local cells for GCRM mesh.");
-  }
-  else {
-    rval = create_local_cells(vertices_on_local_cells, start_vertex, faces);
-    ERRORR(rval, "Failed to create local cells for GCRM mesh.");
-  }
-
-  // Set tag for numCellGroups
-  Tag numCellGroupsTag = 0;
-  rval = mbImpl->tag_get_handle("__NUM_CELL_GROUPS", 1, MB_TYPE_INTEGER, numCellGroupsTag, MB_TAG_SPARSE | MB_TAG_CREAT);
-  ERRORR(rval, "Failed to get __NUM_CELL_GROUPS tag.");
-  rval = mbImpl->tag_set_data(numCellGroupsTag, &_fileSet, 1, &numCellGroups);
-  ERRORR(rval, "Failed to set data for __NUM_CELL_GROUPS tag.");
+  // Create local cells with padding
+  rval = create_padded_local_cells(vertices_on_local_cells, start_vertex, faces);
+  ERRORR(rval, "Failed to create local cells for GCRM mesh.");
 
   if (createGatherSet) {
     EntityHandle gather_set;
@@ -408,15 +389,9 @@ ErrorCode NCHelperGCRM::create_mesh(Range& faces)
       ERRORR(rval, "Failed to create gather set edges for GCRM mesh.");
     }
 
-    // Create gather set cells, either unpadded or padded
-    if (noMixedElements) {
-      rval = create_padded_gather_set_cells(gather_set, start_gather_set_vertex);
-      ERRORR(rval, "Failed to create padded gather set cells for GCRM mesh.");
-    }
-    else {
-      rval = create_gather_set_cells(gather_set, start_gather_set_vertex);
-      ERRORR(rval, "Failed to create gather set cells for GCRM mesh.");
-    }
+    // Create gather set cells with padding
+    rval = create_padded_gather_set_cells(gather_set, start_gather_set_vertex);
+    ERRORR(rval, "Failed to create gather set cells for GCRM mesh.");
   }
 
   return MB_SUCCESS;
@@ -541,19 +516,13 @@ ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset_allocate(std::vector<ReadNC
       }
 
       // Get ptr to tag space
-      if (vdatas[i].entLoc == ReadNC::ENTLOCFACE && numCellGroups > 1) {
-        // For a cell variable that is NOT on one contiguous chunk of faces, defer its tag space allocation
-        vdatas[i].varDatas[t] = NULL;
-      }
-      else {
-        assert(1 == range->psize());
-        void* data;
-        int count;
-        rval = mbImpl->tag_iterate(vdatas[i].varTags[t], range->begin(), range->end(), count, data);
-        ERRORR(rval, "Failed to iterate tag.");
-        assert((unsigned)count == range->size());
-        vdatas[i].varDatas[t] = data;
-      }
+      assert(1 == range->psize());
+      void* data;
+      int count;
+      rval = mbImpl->tag_iterate(vdatas[i].varTags[t], range->begin(), range->end(), count, data);
+      ERRORR(rval, "Failed to iterate tag.");
+      assert((unsigned)count == range->size());
+      vdatas[i].varDatas[t] = data;
     }
   }
 
@@ -563,7 +532,6 @@ ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset_allocate(std::vector<ReadNC
 #ifdef PNETCDF_FILE
 ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset_async(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums)
 {
-  Interface*& mbImpl = _readNC->mbImpl;
   bool& noEdges = _readNC->noEdges;
   DebugOutput& dbgOut = _readNC->dbgOut;
 
@@ -647,32 +615,9 @@ ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset_async(std::vector<ReadNC::V
           success = NCFUNC(wait_all)(_fileId, requests.size(), &requests[0], &statuss[0]);
           ERRORS(success, "Failed on wait_all.");
 
-          if (vdatas[i].entLoc == ReadNC::ENTLOCFACE && numCellGroups > 1) {
-            // For a cell variable that is NOT on one contiguous chunk of faces, allocate tag space for
-            // each cell group, and utilize cellHandleToGlobalID map to read tag data
-            Range::iterator iter = facesOwned.begin();
-            while (iter != facesOwned.end()) {
-              int count;
-              void* ptr;
-              rval = mbImpl->tag_iterate(vdatas[i].varTags[t], iter, facesOwned.end(), count, ptr);
-              ERRORR(rval, "Failed to iterate tag on owned faces.");
-
-              for (int j = 0; j < count; j++) {
-                int global_cell_idx = cellHandleToGlobalID[*(iter + j)]; // Global cell index, 1 based
-                int local_cell_idx = localGidCells.index(global_cell_idx); // Local cell index, 0 based
-                assert(local_cell_idx != -1);
-                for (int level = 0; level < vdatas[i].numLev; level++)
-                  ((double*) ptr)[j * vdatas[i].numLev + level] = tmpdoubledata[local_cell_idx * vdatas[i].numLev + level];
-              }
-
-              iter += count;
-            }
-          }
-          else {
-            void* data = vdatas[i].varDatas[t];
-            for (std::size_t idx = 0; idx != tmpdoubledata.size(); idx++)
-              ((double*) data)[idx] = tmpdoubledata[idx];
-          }
+          void* data = vdatas[i].varDatas[t];
+          for (std::size_t idx = 0; idx != tmpdoubledata.size(); idx++)
+            ((double*) data)[idx] = tmpdoubledata[idx];
 
           break;
         }
@@ -694,17 +639,6 @@ ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset_async(std::vector<ReadNC::V
     }
   }
 
-  for (unsigned int i = 0; i < vdatas.size(); i++) {
-    if (noEdges && vdatas[i].entLoc == ReadNC::ENTLOCEDGE)
-      continue;
-
-    /*for (unsigned int t = 0; t < tstep_nums.size(); t++) {
-      dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
-      ErrorCode tmp_rval = convert_variable(vdatas[i], t);
-      if (MB_SUCCESS != tmp_rval)
-        rval = tmp_rval;
-    }*/
-  }
   // Debug output, if requested
   if (1 == dbgOut.get_verbosity()) {
     dbgOut.printf(1, "Read variables: %s", vdatas.begin()->varName.c_str());
@@ -790,32 +724,9 @@ ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset(std::vector<ReadNC::VarData
           }
           assert(ic == pLocalGid->psize());
 
-          if (vdatas[i].entLoc == ReadNC::ENTLOCFACE && numCellGroups > 1) {
-            // For a cell variable that is NOT on one contiguous chunk of faces, allocate tag space for
-            // each cell group, and utilize cellHandleToGlobalID map to read tag data
-            Range::iterator iter = facesOwned.begin();
-            while (iter != facesOwned.end()) {
-              int count;
-              void* ptr;
-              rval = mbImpl->tag_iterate(vdatas[i].varTags[t], iter, facesOwned.end(), count, ptr);
-              ERRORR(rval, "Failed to iterate tag on owned faces.");
-
-              for (int j = 0; j < count; j++) {
-                int global_cell_idx = cellHandleToGlobalID[*(iter + j)]; // Global cell index, 1 based
-                int local_cell_idx = localGidCells.index(global_cell_idx); // Local cell index, 0 based
-                assert(local_cell_idx != -1);
-                for (int level = 0; level < vdatas[i].numLev; level++)
-                  ((double*) ptr)[j * vdatas[i].numLev + level] = tmpdoubledata[local_cell_idx * vdatas[i].numLev + level];
-              }
-
-              iter += count;
-            }
-          }
-          else {
-            void* data = vdatas[i].varDatas[t];
-            for (std::size_t idx = 0; idx != tmpdoubledata.size(); idx++)
-              ((double*) data)[idx] = tmpdoubledata[idx];
-          }
+          void* data = vdatas[i].varDatas[t];
+          for (std::size_t idx = 0; idx != tmpdoubledata.size(); idx++)
+            ((double*) data)[idx] = tmpdoubledata[idx];
 
           break;
         }
@@ -836,18 +747,6 @@ ErrorCode NCHelperGCRM::read_ucd_variables_to_nonset(std::vector<ReadNC::VarData
     }
   }
 
-  for (unsigned int i = 0; i < vdatas.size(); i++) {
-    if (noEdges && vdatas[i].entLoc == ReadNC::ENTLOCEDGE)
-      continue;
-
-   /* for (unsigned int t = 0; t < tstep_nums.size(); t++) {
-      dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
-      ErrorCode tmp_rval = convert_variable(vdatas[i], t);
-      if (MB_SUCCESS != tmp_rval)
-        rval = tmp_rval;
-    }*/
-  }
-
   // Debug output, if requested
   if (1 == dbgOut.get_verbosity()) {
     dbgOut.printf(1, "Read variables: %s", vdatas.begin()->varName.c_str());
@@ -865,8 +764,8 @@ ErrorCode NCHelperGCRM::redistribute_local_cells(int start_cell_idx)
   // If possible, apply Zoltan partition
   if (_readNC->partMethod == ScdParData::RCBZOLTAN) {
 #if defined(USE_MPI) && defined(HAVE_ZOLTAN)
-    // Read lat/lon coordinates of cell centers
-    // then convert to spherical , and use them as input to zoltan partition
+    // Read lat/lon coordinates of cell centers, then convert spherical to
+    // Cartesian, and use them as input to Zoltan partition
     int xCellVarId;
     int success = NCFUNC(inq_varid)(_fileId, "grid_center_lat", &xCellVarId);
     ERRORS(success, "Failed to get variable id of grid_center_lat.");
@@ -874,7 +773,7 @@ ErrorCode NCHelperGCRM::redistribute_local_cells(int start_cell_idx)
     NCDF_SIZE read_start = static_cast<NCDF_SIZE>(start_cell_idx - 1);
     NCDF_SIZE read_count = static_cast<NCDF_SIZE>(nLocalCells);
     success = NCFUNCAG(_vara_double)(_fileId, xCellVarId, &read_start, &read_count, &xCell[0]);
-    ERRORS(success, "Failed to read xCell data.");
+    ERRORS(success, "Failed to read grid_center_lat data.");
 
     // Read y coordinates of cell centers
     int yCellVarId;
@@ -882,14 +781,12 @@ ErrorCode NCHelperGCRM::redistribute_local_cells(int start_cell_idx)
     ERRORS(success, "Failed to get variable id of grid_center_lon.");
     std::vector<double> yCell(nLocalCells);
     success = NCFUNCAG(_vara_double)(_fileId, yCellVarId, &read_start, &read_count, &yCell[0]);
-    ERRORS(success, "Failed to read yCell data.");
+    ERRORS(success, "Failed to read grid_center_lon data.");
 
+    // Convert lon/lat/rad to x/y/z
     std::vector<double> zCell(nLocalCells);
-    // convert to xyz cartesian coordinates
-
-    double rad=8000; // this is just approx x is lat, y is lon
-    for (int i=0; i<nLocalCells; i++)
-    {
+    double rad = 8000.0; // This is just a approximate radius
+    for (int i = 0; i < nLocalCells; i++) {
       double cosphi = cos(xCell[i]);
       double zmult = sin(xCell[i]);
       double xmult = cosphi * cos(yCell[i]);
@@ -898,6 +795,7 @@ ErrorCode NCHelperGCRM::redistribute_local_cells(int start_cell_idx)
       yCell[i] = rad * ymult;
       zCell[i] = rad * zmult;
     }
+
     // Zoltan partition using RCB; maybe more studies would be good, as to which partition
     // is better
     Interface*& mbImpl = _readNC->mbImpl;
@@ -991,20 +889,18 @@ ErrorCode NCHelperGCRM::create_local_vertices(const std::vector<int>& vertices_o
     ERRORR(MB_FAILURE, "Couldn't find 'layers' or 'interfaces' variable.");
   }
 
-  {
-    // Decide whether down is positive
-    char posval[10] = {0};
-    int success = NCFUNC(get_att_text)(_fileId, (*vmit).second.varId, "positive", posval);
-    if (0 == success && !strncmp(posval, "down", 4)) {
-      for (std::vector<double>::iterator dvit = levVals.begin(); dvit != levVals.end(); ++dvit)
-        (*dvit) *= -1.0;
-    }
+  // Decide whether down is positive
+  char posval[10] = {0};
+  int success = NCFUNC(get_att_text)(_fileId, (*vmit).second.varId, "positive", posval);
+  if (0 == success && !strncmp(posval, "down", 4)) {
+    for (std::vector<double>::iterator dvit = levVals.begin(); dvit != levVals.end(); ++dvit)
+      (*dvit) *= -1.0;
   }
 
   // Read x coordinates for local vertices
   double* xptr = arrays[0];
   int xVertexVarId;
-  int success = NCFUNC(inq_varid)(_fileId, "grid_corner_lon", &xVertexVarId);
+  success = NCFUNC(inq_varid)(_fileId, "grid_corner_lon", &xVertexVarId);
   ERRORS(success, "Failed to get variable id of grid_corner_lon.");
   size_t indexInArray = 0;
   for (Range::pair_iterator pair_iter = localGidVerts.pair_begin();
@@ -1074,11 +970,10 @@ ErrorCode NCHelperGCRM::create_local_vertices(const std::vector<int>& vertices_o
 
   // Convert lon/lat/rad to x/y/z
   double* zptr = arrays[2];
-  //const double pideg = acos(-1.0) / 180.0;
   double rad = 8000.0 + levVals[0];
   for (int i = 0; i < nLocalVertices; i++) {
     double cosphi = cos(yptr[i]);
-    double zmult =  sin(yptr[i]);
+    double zmult = sin(yptr[i]);
     double xmult = cosphi * cos(xptr[i]);
     double ymult = cosphi * sin(xptr[i]);
     xptr[i] = rad * xmult;
@@ -1100,7 +995,7 @@ ErrorCode NCHelperGCRM::create_local_edges(EntityHandle start_vertex)
   int success = NCFUNC(inq_varid)(_fileId, "cell_edges", &edgesOnCellVarId);
   ERRORS(success, "Failed to get variable id of cell_edges.");
 
-  std::vector<int> edges_on_local_cells(nLocalCells * maxEdgesPerCell);
+  std::vector<int> edges_on_local_cells(nLocalCells * EDGES_PER_CELL);
   dbgOut.tprintf(1, "   edges_on_local_cells.size() = %d\n", (int)edges_on_local_cells.size());
 
 #ifdef PNETCDF_FILE
@@ -1118,7 +1013,7 @@ ErrorCode NCHelperGCRM::create_local_edges(EntityHandle start_vertex)
     dbgOut.tprintf(1, "   starth = %d\n", (int)starth);
     dbgOut.tprintf(1, "   endh = %d\n", (int)endh);
     NCDF_SIZE read_starts[2] = {static_cast<NCDF_SIZE>(starth - 1), 0};
-    NCDF_SIZE read_counts[2] = {static_cast<NCDF_SIZE>(endh - starth + 1), static_cast<NCDF_SIZE>(maxEdgesPerCell)};
+    NCDF_SIZE read_counts[2] = {static_cast<NCDF_SIZE>(endh - starth + 1), static_cast<NCDF_SIZE>(EDGES_PER_CELL)};
 
     // Do a partial read in each subrange
 #ifdef PNETCDF_FILE
@@ -1131,7 +1026,7 @@ ErrorCode NCHelperGCRM::create_local_edges(EntityHandle start_vertex)
     ERRORS(success, "Failed to read cell_edges data in a loop");
 
     // Increment the index for next subrange
-    indexInArray += (endh - starth + 1) * maxEdgesPerCell;
+    indexInArray += (endh - starth + 1) * EDGES_PER_CELL;
   }
 
 #ifdef PNETCDF_FILE
@@ -1141,9 +1036,8 @@ ErrorCode NCHelperGCRM::create_local_edges(EntityHandle start_vertex)
 #endif
 
   // GCRM is 0 based, convert edge indices from 0 to 1 based
-  for (std::size_t idx = 0; idx < edges_on_local_cells.size(); idx++) {
-      edges_on_local_cells[idx] += 1;
-  }
+  for (std::size_t idx = 0; idx < edges_on_local_cells.size(); idx++)
+    edges_on_local_cells[idx] += 1;
 
   // Collect local edges
   std::sort(edges_on_local_cells.begin(), edges_on_local_cells.end());
@@ -1217,16 +1111,12 @@ ErrorCode NCHelperGCRM::create_local_edges(EntityHandle start_vertex)
   ERRORS(success, "Failed on wait_all.");
 #endif
 
-  // GCRM is 0 based, convert edge indices from 0 to 1 based
-  for (int idx = 0; idx < nLocalEdges*2; idx++) {
-      vertices_on_local_edges[idx] += 1;
-  }
-
   // Populate connectivity data for local edges
   // Convert in-place from int (stored in the first half) to EntityHandle
   // Reading backward is the trick
   for (int edge_vert = nLocalEdges * 2 - 1; edge_vert >= 0; edge_vert--) {
-    int global_vert_idx = vertices_on_local_edges[edge_vert]; // Global vertex index, 1 based
+    // Note, indices stored in vertices_on_local_edges are 0 based
+    int global_vert_idx = vertices_on_local_edges[edge_vert] + 1; // Global vertex index, 1 based
     int local_vert_idx = localGidVerts.index(global_vert_idx); // Local vertex index, 0 based
     assert(local_vert_idx != -1);
     conn_arr_edges[edge_vert] = start_vertex + local_vert_idx;
@@ -1235,104 +1125,16 @@ ErrorCode NCHelperGCRM::create_local_edges(EntityHandle start_vertex)
   return MB_SUCCESS;
 }
 
-ErrorCode NCHelperGCRM::create_local_cells(const std::vector<int>& vertices_on_local_cells,
-                                                    EntityHandle start_vertex, Range& faces)
-{
-  Interface*& mbImpl = _readNC->mbImpl;
-  Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
-
-  // Divide local cells into groups based on the number of edges
-  Range local_cells_with_n_edges[DEFAULT_MAX_EDGES_PER_CELL + 1];
-  // Insert larger values before smaller ones to increase efficiency
-  for (int i = nLocalCells - 1; i >= 0; i--) {
-    int num_edges = DEFAULT_MAX_EDGES_PER_CELL;
-    local_cells_with_n_edges[num_edges].insert(localGidCells[i]); // Global cell index
-  }
-
-  std::vector<int> num_edges_on_cell_groups;
-  for (int i = 3; i <= maxEdgesPerCell; i++) {
-    if (local_cells_with_n_edges[i].size() > 0)
-      num_edges_on_cell_groups.push_back(i);
-  }
-  numCellGroups = num_edges_on_cell_groups.size();
-
-  EntityHandle* conn_arr_local_cells_with_n_edges[DEFAULT_MAX_EDGES_PER_CELL + 1];
-  for (int i = 0; i < numCellGroups; i++) {
-    int num_edges_per_cell = num_edges_on_cell_groups[i];
-    int num_group_cells = (int)local_cells_with_n_edges[num_edges_per_cell].size();
-
-    // Create local cells for each non-empty cell group
-    EntityHandle start_element;
-    ErrorCode rval = _readNC->readMeshIface->get_element_connect(num_group_cells, num_edges_per_cell, MBPOLYGON, 0, start_element,
-                                                       conn_arr_local_cells_with_n_edges[num_edges_per_cell], num_group_cells);
-    ERRORR(rval, "Failed to create cells");
-    faces.insert(start_element, start_element + num_group_cells - 1);
-
-    // Add local cells to the file set
-    Range local_cells_range(start_element, start_element + num_group_cells - 1);
-    rval = _readNC->mbImpl->add_entities(_fileSet, local_cells_range);
-    ERRORR(rval, "Failed to add local cells to the file set.");
-
-    // Get ptr to gid memory for local cells
-    int count = 0;
-    void* data = NULL;
-    rval = mbImpl->tag_iterate(mGlobalIdTag, local_cells_range.begin(), local_cells_range.end(), count, data);
-    ERRORR(rval, "Failed to iterate global id tag on local cells.");
-    assert(count == num_group_cells);
-    int* gid_data = (int*) data;
-    std::copy(local_cells_with_n_edges[num_edges_per_cell].begin(), local_cells_with_n_edges[num_edges_per_cell].end(), gid_data);
-
-    // Set connectivity array with proper local vertices handles
-    for (int j = 0; j < num_group_cells; j++) {
-      EntityHandle global_cell_idx = local_cells_with_n_edges[num_edges_per_cell][j]; // Global cell index, 1 based
-      int local_cell_idx = localGidCells.index(global_cell_idx); // Local cell index, 0 based
-      assert(local_cell_idx != -1);
-
-      if (numCellGroups > 1) {
-        // Populate cellHandleToGlobalID map to read cell variables
-        cellHandleToGlobalID[start_element + j] = global_cell_idx;
-      }
-
-      for (int k = 0; k < num_edges_per_cell; k++) {
-        EntityHandle global_vert_idx = vertices_on_local_cells[local_cell_idx * maxEdgesPerCell + k]; // Global vertex index, 1 based
-        int local_vert_idx = localGidVerts.index(global_vert_idx); // Local vertex index, 0 based
-        assert(local_vert_idx != -1);
-        conn_arr_local_cells_with_n_edges[num_edges_per_cell][j * num_edges_per_cell + k] =
-            start_vertex + local_vert_idx;
-      }
-      // make sure that if some nodes are repeated, they are at the end of the connectivity array
-      // so, pentagons as hexagons should have a connectivity like 123455 and not 122345
-      EntityHandle *pvertex= &(conn_arr_local_cells_with_n_edges[num_edges_per_cell][j * num_edges_per_cell ]);
-      for (int  k = 0; k < num_edges_per_cell-2; k++)
-      {
-        if( *(pvertex+k) == *(pvertex+k+1) )
-        {
-          // shift the connectivity
-          for (int kk=k+1; kk<num_edges_per_cell-1; kk++)
-          {
-            *(pvertex+kk)=*(pvertex+kk+1);
-          }
-        }
-      }
-    }
-  }
-
-  return MB_SUCCESS;
-}
-
 ErrorCode NCHelperGCRM::create_padded_local_cells(const std::vector<int>& vertices_on_local_cells,
                                                   EntityHandle start_vertex, Range& faces)
 {
   Interface*& mbImpl = _readNC->mbImpl;
   Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
 
-  // Only one group of cells (each cell is represented by a polygon with maxEdgesPerCell edges)
-  numCellGroups = 1;
-
-  // Create cells for this cell group
+  // Create cells
   EntityHandle start_element;
   EntityHandle* conn_arr_local_cells = NULL;
-  ErrorCode rval = _readNC->readMeshIface->get_element_connect(nLocalCells, maxEdgesPerCell, MBPOLYGON, 0, start_element, conn_arr_local_cells,
+  ErrorCode rval = _readNC->readMeshIface->get_element_connect(nLocalCells, EDGES_PER_CELL, MBPOLYGON, 0, start_element, conn_arr_local_cells,
                                                     // Might have to create gather mesh later
                                                     (createGatherSet ? nLocalCells + nCells : nLocalCells));
   ERRORR(rval, "Failed to create cells.");
@@ -1353,14 +1155,15 @@ ErrorCode NCHelperGCRM::create_padded_local_cells(const std::vector<int>& vertic
   std::copy(localGidCells.begin(), localGidCells.end(), gid_data);
 
   // Set connectivity array with proper local vertices handles
-  // vertices_on_local_cells array was already corrected to have the last vertices padded
-  // no need for extra checks considering
+  // vertices_on_local_cells array was already corrected to have
+  // the last vertices repeated for pentagons, e.g. 122345 => 123455
   for (int local_cell_idx = 0; local_cell_idx < nLocalCells; local_cell_idx++) {
-    for (int i = 0; i < maxEdgesPerCell; i++) {
-      EntityHandle global_vert_idx = vertices_on_local_cells[local_cell_idx * maxEdgesPerCell + i]; // Global vertex index, 1 based
+    for (int i = 0; i < EDGES_PER_CELL; i++) {
+      // Note, indices stored in vertices_on_local_cells are 1 based
+      EntityHandle global_vert_idx = vertices_on_local_cells[local_cell_idx * EDGES_PER_CELL + i]; // Global vertex index, 1 based
       int local_vert_idx = localGidVerts.index(global_vert_idx); // Local vertex index, 0 based
       assert(local_vert_idx != -1);
-      conn_arr_local_cells[local_cell_idx * maxEdgesPerCell + i] = start_vertex + local_vert_idx;
+      conn_arr_local_cells[local_cell_idx * EDGES_PER_CELL + i] = start_vertex + local_vert_idx;
     }
   }
 
@@ -1387,8 +1190,8 @@ ErrorCode NCHelperGCRM::create_gather_set_vertices(EntityHandle gather_set, Enti
   // Read x coordinates for gather set vertices
   double* xptr = arrays[0];
   int xVertexVarId;
-  int success = NCFUNC(inq_varid)(_fileId, "xVertex", &xVertexVarId);
-  ERRORS(success, "Failed to get variable id of xVertex.");
+  int success = NCFUNC(inq_varid)(_fileId, "grid_corner_lon", &xVertexVarId);
+  ERRORS(success, "Failed to get variable id of grid_corner_lon.");
   NCDF_SIZE read_start = 0;
   NCDF_SIZE read_count = static_cast<NCDF_SIZE>(nVertices);
 #ifdef PNETCDF_FILE
@@ -1396,49 +1199,44 @@ ErrorCode NCHelperGCRM::create_gather_set_vertices(EntityHandle gather_set, Enti
   success = NCFUNC(begin_indep_data)(_fileId);
   ERRORS(success, "Failed to begin independent I/O mode.");
   success = NCFUNCG(_vara_double)(_fileId, xVertexVarId, &read_start, &read_count, xptr);
-  ERRORS(success, "Failed to read xVertex data.");
+  ERRORS(success, "Failed to read grid_corner_lon data.");
   success = NCFUNC(end_indep_data)(_fileId);
   ERRORS(success, "Failed to end independent I/O mode.");
 #else
   success = NCFUNCG(_vara_double)(_fileId, xVertexVarId, &read_start, &read_count, xptr);
-  ERRORS(success, "Failed to read xVertex data.");
+  ERRORS(success, "Failed to read grid_corner_lon data.");
 #endif
 
   // Read y coordinates for gather set vertices
   double* yptr = arrays[1];
   int yVertexVarId;
-  success = NCFUNC(inq_varid)(_fileId, "yVertex", &yVertexVarId);
-  ERRORS(success, "Failed to get variable id of yVertex.");
+  success = NCFUNC(inq_varid)(_fileId, "grid_corner_lat", &yVertexVarId);
+  ERRORS(success, "Failed to get variable id of grid_corner_lat.");
 #ifdef PNETCDF_FILE
   // Enter independent I/O mode, since this read is only for the gather processor
   success = NCFUNC(begin_indep_data)(_fileId);
   ERRORS(success, "Failed to begin independent I/O mode.");
   success = NCFUNCG(_vara_double)(_fileId, yVertexVarId, &read_start, &read_count, yptr);
-  ERRORS(success, "Failed to read yVertex data.");
+  ERRORS(success, "Failed to read grid_corner_lat data.");
   success = NCFUNC(end_indep_data)(_fileId);
   ERRORS(success, "Failed to end independent I/O mode.");
 #else
   success = NCFUNCG(_vara_double)(_fileId, yVertexVarId, &read_start, &read_count, yptr);
-  ERRORS(success, "Failed to read yVertex data.");
+  ERRORS(success, "Failed to read grid_corner_lat data.");
 #endif
 
-  // Read z coordinates for gather set vertices
+  // Convert lon/lat/rad to x/y/z
   double* zptr = arrays[2];
-  int zVertexVarId;
-  success = NCFUNC(inq_varid)(_fileId, "zVertex", &zVertexVarId);
-  ERRORS(success, "Failed to get variable id of zVertex.");
-#ifdef PNETCDF_FILE
-  // Enter independent I/O mode, since this read is only for the gather processor
-  success = NCFUNC(begin_indep_data)(_fileId);
-  ERRORS(success, "Failed to begin independent I/O mode.");
-  success = NCFUNCG(_vara_double)(_fileId, zVertexVarId, &read_start, &read_count, zptr);
-  ERRORS(success, "Failed to read zVertex data.");
-  success = NCFUNC(end_indep_data)(_fileId);
-  ERRORS(success, "Failed to end independent I/O mode.");
-#else
-  success = NCFUNCG(_vara_double)(_fileId, zVertexVarId, &read_start, &read_count, zptr);
-  ERRORS(success, "Failed to read zVertex data.");
-#endif
+  double rad = 8000.0 + levVals[0];
+  for (int i = 0; i < nVertices; i++) {
+    double cosphi = cos(yptr[i]);
+    double zmult = sin(yptr[i]);
+    double xmult = cosphi * cos(xptr[i]);
+    double ymult = cosphi * sin(xptr[i]);
+    xptr[i] = rad * xmult;
+    yptr[i] = rad * ymult;
+    zptr[i] = rad * zmult;
+  }
 
   // Get ptr to GID memory for gather set vertices
   int count = 0;
@@ -1504,8 +1302,8 @@ ErrorCode NCHelperGCRM::create_gather_set_edges(EntityHandle gather_set, EntityH
    // Convert in-place from int (stored in the first half) to EntityHandle
    // Reading backward is the trick
    for (int edge_vert = nEdges * 2 - 1; edge_vert >= 0; edge_vert--) {
-     int gather_set_vert_idx = vertices_on_gather_set_edges[edge_vert]; // Global vertex index, 1 based
-     gather_set_vert_idx--; // 1 based -> 0 based
+     // Note, indices stored in vertices_on_gather_set_edges are 0 based
+     int gather_set_vert_idx = vertices_on_gather_set_edges[edge_vert]; // Global vertex index, 0 based
      // Connectivity array is shifted by where the gather set vertices start
      conn_arr_gather_set_edges[edge_vert] = gather_set_start_vertex + gather_set_vert_idx;
    }
@@ -1513,121 +1311,15 @@ ErrorCode NCHelperGCRM::create_gather_set_edges(EntityHandle gather_set, EntityH
    return MB_SUCCESS;
 }
 
-ErrorCode NCHelperGCRM::create_gather_set_cells(EntityHandle gather_set, EntityHandle gather_set_start_vertex)
-{
-  Interface*& mbImpl = _readNC->mbImpl;
-
-  // Read number of edges on each gather set cell
-  int nEdgesOnCellVarId;
-  int success = NCFUNC(inq_varid)(_fileId, "nEdgesOnCell", &nEdgesOnCellVarId);
-  ERRORS(success, "Failed to get variable id of nEdgesOnCell.");
-  std::vector<int> num_edges_on_gather_set_cells(nCells);
-  NCDF_SIZE read_start = 0;
-  NCDF_SIZE read_count = static_cast<NCDF_SIZE>(nCells);
-#ifdef PNETCDF_FILE
-  // Enter independent I/O mode, since this read is only for the gather processor
-  success = NCFUNC(begin_indep_data)(_fileId);
-  ERRORS(success, "Failed to begin independent I/O mode.");
-  success = NCFUNCG(_vara_int)(_fileId, nEdgesOnCellVarId, &read_start, &read_count, &num_edges_on_gather_set_cells[0]);
-  ERRORS(success, "Failed to read nEdgesOnCell data.");
-  success = NCFUNC(end_indep_data)(_fileId);
-  ERRORS(success, "Failed to end independent I/O mode.");
-#else
-  success = NCFUNCG(_vara_int)(_fileId, nEdgesOnCellVarId, &read_start, &read_count, &num_edges_on_gather_set_cells[0]);
-  ERRORS(success, "Failed to read nEdgesOnCell data.");
-#endif
-
-  // Read vertices on each gather set cell (connectivity)
-  int verticesOnCellVarId;
-  success = NCFUNC(inq_varid)(_fileId, "verticesOnCell", &verticesOnCellVarId);
-  ERRORS(success, "Failed to get variable id of verticesOnCell.");
-  std::vector<int> vertices_on_gather_set_cells(nCells * maxEdgesPerCell);
-  NCDF_SIZE read_starts[2] = {0, 0};
-  NCDF_SIZE read_counts[2] = {static_cast<NCDF_SIZE>(nCells), static_cast<NCDF_SIZE>(maxEdgesPerCell)};
-#ifdef PNETCDF_FILE
-  // Enter independent I/O mode, since this read is only for the gather processor
-  success = NCFUNC(begin_indep_data)(_fileId);
-  ERRORS(success, "Failed to begin independent I/O mode.");
-  success = NCFUNCG(_vara_int)(_fileId, verticesOnCellVarId, read_starts, read_counts, &vertices_on_gather_set_cells[0]);
-  ERRORS(success, "Failed to read verticesOnCell data.");
-  success = NCFUNC(end_indep_data)(_fileId);
-  ERRORS(success, "Failed to end independent I/O mode.");
-#else
-  success = NCFUNCG(_vara_int)(_fileId, verticesOnCellVarId, read_starts, read_counts, &vertices_on_gather_set_cells[0]);
-  ERRORS(success, "Failed to read verticesOnCell data.");
-#endif
-
-  // Divide gather set cells into groups based on the number of edges
-  Range gather_set_cells_with_n_edges[DEFAULT_MAX_EDGES_PER_CELL + 1];
-  // Insert larger values before smaller values to increase efficiency
-  for (int i = nCells - 1; i >= 0; i--) {
-    int num_edges = num_edges_on_gather_set_cells[i];
-    gather_set_cells_with_n_edges[num_edges].insert(i + 1); // 0 based -> 1 based
-  }
-
-  // Create gather set cells
-  EntityHandle* conn_arr_gather_set_cells_with_n_edges[DEFAULT_MAX_EDGES_PER_CELL + 1];
-  for (int num_edges_per_cell = 3; num_edges_per_cell <= maxEdgesPerCell; num_edges_per_cell++) {
-    int num_group_cells = gather_set_cells_with_n_edges[num_edges_per_cell].size();
-    if (num_group_cells > 0) {
-      EntityHandle start_element;
-      ErrorCode rval = _readNC->readMeshIface->get_element_connect(num_group_cells, num_edges_per_cell, MBPOLYGON, 0, start_element,
-                                                         conn_arr_gather_set_cells_with_n_edges[num_edges_per_cell], num_group_cells);
-      ERRORR(rval, "Failed to create cells.");
-
-      // Add cells to the gather set
-      Range gather_set_cells_range(start_element, start_element + num_group_cells - 1);
-      rval = mbImpl->add_entities(gather_set, gather_set_cells_range);
-      ERRORR(rval, "Failed to add cells to the gather set.");
-
-      for (int j = 0; j < num_group_cells; j++) {
-        int gather_set_cell_idx = gather_set_cells_with_n_edges[num_edges_per_cell][j]; // Global cell index, 1 based
-        gather_set_cell_idx--; // 1 based -> 0 based
-
-        for (int k = 0; k < num_edges_per_cell; k++) {
-          EntityHandle gather_set_vert_idx = vertices_on_gather_set_cells[gather_set_cell_idx * maxEdgesPerCell + k]; // Global vertex index, 1 based
-          gather_set_vert_idx--; // 1 based -> 0 based
-
-          // Connectivity array is shifted by where the gather set vertices start
-          conn_arr_gather_set_cells_with_n_edges[num_edges_per_cell][j * num_edges_per_cell + k] =
-            gather_set_start_vertex + gather_set_vert_idx;
-        }
-      }
-    }
-  }
-
-  return MB_SUCCESS;
-}
-
 ErrorCode NCHelperGCRM::create_padded_gather_set_cells(EntityHandle gather_set, EntityHandle gather_set_start_vertex)
 {
   Interface*& mbImpl = _readNC->mbImpl;
 
-  // Read number of edges on each gather set cell
-  int nEdgesOnCellVarId;
-  int success = NCFUNC(inq_varid)(_fileId, "nEdgesOnCell", &nEdgesOnCellVarId);
-  ERRORS(success, "Failed to get variable id of nEdgesOnCell.");
-  std::vector<int> num_edges_on_gather_set_cells(nCells);
-  NCDF_SIZE read_start = 0;
-  NCDF_SIZE read_count = static_cast<NCDF_SIZE>(nCells);
-#ifdef PNETCDF_FILE
-  // Enter independent I/O mode, since this read is only for the gather processor
-  success = NCFUNC(begin_indep_data)(_fileId);
-  ERRORS(success, "Failed to begin independent I/O mode.");
-  success = NCFUNCG(_vara_int)(_fileId, nEdgesOnCellVarId, &read_start, &read_count, &num_edges_on_gather_set_cells[0]);
-  ERRORS(success, "Failed to read nEdgesOnCell data.");
-  success = NCFUNC(end_indep_data)(_fileId);
-  ERRORS(success, "Failed to end independent I/O mode.");
-#else
-  success = NCFUNCG(_vara_int)(_fileId, nEdgesOnCellVarId, &read_start, &read_count, &num_edges_on_gather_set_cells[0]);
-  ERRORS(success, "Failed to read nEdgesOnCell data.");
-#endif
-
   // Create gather set cells
   EntityHandle start_element;
   EntityHandle* conn_arr_gather_set_cells = NULL;
   // Don't need to specify allocation number here, because we know enough cells were created before
-  ErrorCode rval = _readNC->readMeshIface->get_element_connect(nCells, maxEdgesPerCell, MBPOLYGON, 0, start_element, conn_arr_gather_set_cells);
+  ErrorCode rval = _readNC->readMeshIface->get_element_connect(nCells, EDGES_PER_CELL, MBPOLYGON, 0, start_element, conn_arr_gather_set_cells);
   ERRORR(rval, "Failed to create cells.");
 
   // Add cells to the gather set
@@ -1637,41 +1329,46 @@ ErrorCode NCHelperGCRM::create_padded_gather_set_cells(EntityHandle gather_set,
 
   // Read vertices on each gather set cell (connectivity)
   int verticesOnCellVarId;
-  success = NCFUNC(inq_varid)(_fileId, "verticesOnCell", &verticesOnCellVarId);
-  ERRORS(success, "Failed to get variable id of verticesOnCell.");
+  int success = NCFUNC(inq_varid)(_fileId, "cell_corners", &verticesOnCellVarId);
+  ERRORS(success, "Failed to get variable id of cell_corners.");
   // Utilize the memory storage pointed by conn_arr_gather_set_cells
   int* vertices_on_gather_set_cells = (int*) conn_arr_gather_set_cells;
   NCDF_SIZE read_starts[2] = {0, 0};
-  NCDF_SIZE read_counts[2] = {static_cast<NCDF_SIZE>(nCells), static_cast<NCDF_SIZE>(maxEdgesPerCell)};
+  NCDF_SIZE read_counts[2] = {static_cast<NCDF_SIZE>(nCells), static_cast<NCDF_SIZE>(EDGES_PER_CELL)};
 #ifdef PNETCDF_FILE
   // Enter independent I/O mode, since this read is only for the gather processor
   success = NCFUNC(begin_indep_data)(_fileId);
   ERRORS(success, "Failed to begin independent I/O mode.");
   success = NCFUNCG(_vara_int)(_fileId, verticesOnCellVarId, read_starts, read_counts, vertices_on_gather_set_cells);
-  ERRORS(success, "Failed to read verticesOnCell data.");
+  ERRORS(success, "Failed to read cell_corners data.");
   success = NCFUNC(end_indep_data)(_fileId);
   ERRORS(success, "Failed to end independent I/O mode.");
 #else
   success = NCFUNCG(_vara_int)(_fileId, verticesOnCellVarId, read_starts, read_counts, vertices_on_gather_set_cells);
-  ERRORS(success, "Failed to read verticesOnCell data.");
+  ERRORS(success, "Failed to read cell_corners data.");
 #endif
 
-  // Correct gather set cell vertices array in the same way as local cell vertices array,
-  // replace the padded vertices with the last vertices in the corresponding cells
+  // Correct gather set cell vertices array in the same way as local cell vertices array
+  // Pentagons as hexagons should have a connectivity like 123455 and not 122345
   for (int gather_set_cell_idx = 0; gather_set_cell_idx < nCells; gather_set_cell_idx++) {
-    int num_edges = num_edges_on_gather_set_cells[gather_set_cell_idx];
-    int idx_in_gather_set_vert_arr = gather_set_cell_idx * maxEdgesPerCell;
-    int last_vert_idx = vertices_on_gather_set_cells[idx_in_gather_set_vert_arr + num_edges - 1];
-    for (int i = num_edges; i < maxEdgesPerCell; i++)
-      vertices_on_gather_set_cells[idx_in_gather_set_vert_arr + i] = last_vert_idx;
+    int* pvertex = vertices_on_gather_set_cells + gather_set_cell_idx * EDGES_PER_CELL;
+    for (int k = 0; k < EDGES_PER_CELL - 2; k++) {
+      if (*(pvertex + k) == *(pvertex + k + 1)) {
+        // Shift the connectivity
+        for (int kk = k + 1; kk < EDGES_PER_CELL - 1; kk++)
+          *(pvertex + kk) = *(pvertex + kk + 1);
+        // No need to try next k
+        break;
+      }
+    }
   }
 
   // Populate connectivity data for gather set cells
   // Convert in-place from int (stored in the first half) to EntityHandle
   // Reading backward is the trick
-  for (int cell_vert = nCells * maxEdgesPerCell - 1; cell_vert >= 0; cell_vert--) {
-    int gather_set_vert_idx = vertices_on_gather_set_cells[cell_vert]; // Global vertex index, 1 based
-    gather_set_vert_idx--; // 1 based -> 0 based
+  for (int cell_vert = nCells * EDGES_PER_CELL - 1; cell_vert >= 0; cell_vert--) {
+    // Note, indices stored in vertices_on_gather_set_cells are 0 based
+    int gather_set_vert_idx = vertices_on_gather_set_cells[cell_vert]; // Global vertex index, 0 based
     // Connectivity array is shifted by where the gather set vertices start
     conn_arr_gather_set_cells[cell_vert] = gather_set_start_vertex + gather_set_vert_idx;
   }

diff --git a/src/io/NCHelperGCRM.hpp b/src/io/NCHelperGCRM.hpp
index b415bd2..2887e86 100644
--- a/src/io/NCHelperGCRM.hpp
+++ b/src/io/NCHelperGCRM.hpp
@@ -52,11 +52,7 @@ private:
   //! Create local edges (optional)
   ErrorCode create_local_edges(EntityHandle start_vertex);
 
-  //! Create local cells without padding (cells are divided into groups based on the number of edges)
-  ErrorCode create_local_cells(const std::vector<int>& vertices_on_local_cells,
-                                        EntityHandle start_vertex, Range& faces);
-
-  //! Create local cells with padding (padded cells will have the same number of edges)
+  //! Create local cells with padding (pentagons are padded to hexagons)
   ErrorCode create_padded_local_cells(const std::vector<int>& vertices_on_local_cells,
                                       EntityHandle start_vertex, Range& faces);
 
@@ -66,17 +62,11 @@ private:
   //! Create gather set edges (optional)
   ErrorCode create_gather_set_edges(EntityHandle gather_set, EntityHandle gather_set_start_vertex);
 
-  //! Create gather set cells without padding (cells are divided into groups based on the number of edges)
-  ErrorCode create_gather_set_cells(EntityHandle gather_set, EntityHandle gather_set_start_vertex);
-
-  //! Create gather set cells with padding (padded cells will have the same number of edges)
+  //! Create gather set cells with padding (pentagons are padded to hexagons)
   ErrorCode create_padded_gather_set_cells(EntityHandle gather_set, EntityHandle gather_set_start_vertex);
 
 private:
-  int maxEdgesPerCell;
-  int numCellGroups;
   bool createGatherSet;
-  std::map<EntityHandle, int> cellHandleToGlobalID;
   Range facesOwned;
 };
 

diff --git a/test/io/read_gcrm_nc.cpp b/test/io/read_gcrm_nc.cpp
index 35d6bcb..c4a2290 100644
--- a/test/io/read_gcrm_nc.cpp
+++ b/test/io/read_gcrm_nc.cpp
@@ -48,7 +48,7 @@ int main(int argc, char* argv[])
   result += RUN_TEST(test_read_nomesh);
   result += RUN_TEST(test_read_novars);
   result += RUN_TEST(test_read_no_edges);
-  //result += RUN_TEST(test_gather_onevar);
+  result += RUN_TEST(test_gather_onevar);
 
 #ifdef USE_MPI
   fail = MPI_Finalize();
@@ -494,7 +494,80 @@ void test_read_no_edges()
 
 void test_gather_onevar()
 {
-  // TBD
+  Core moab;
+  Interface& mb = moab;
+
+  EntityHandle file_set;
+  ErrorCode rval = mb.create_meshset(MESHSET_SET, file_set);
+  CHECK_ERR(rval);
+
+  std::string opts;
+  get_options(opts);
+
+  // Read cell variable vorticity and create gather set on processor 0
+  opts += ";VARIABLE=vorticity;GATHER_SET=0";
+  rval = mb.load_file(example, &file_set, opts.c_str());
+  CHECK_ERR(rval);
+
+#ifdef USE_MPI
+  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
+  int rank = pcomm->proc_config().proc_rank();
+
+  Range cells, cells_owned;
+  rval = mb.get_entities_by_type(file_set, MBPOLYGON, cells);
+  CHECK_ERR(rval);
+
+  // Get local owned cells
+  rval = pcomm->filter_pstatus(cells, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &cells_owned);
+  CHECK_ERR(rval);
+
+  EntityHandle gather_set = 0;
+  if (0 == rank) {
+    // Get gather set
+    ReadUtilIface* readUtilIface;
+    mb.query_interface(readUtilIface);
+    rval = readUtilIface->get_gather_set(gather_set);
+    CHECK_ERR(rval);
+    assert(gather_set != 0);
+  }
+
+  Tag vorticity_tag0, gid_tag;
+  rval = mb.tag_get_handle("vorticity0", layers, MB_TYPE_DOUBLE, vorticity_tag0, MB_TAG_DENSE);
+  CHECK_ERR(rval);
+
+  rval = mb.tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gid_tag, MB_TAG_DENSE);
+  CHECK_ERR(rval);
+
+  pcomm->gather_data(cells_owned, vorticity_tag0, gid_tag, gather_set, 0);
+
+  if (0 == rank) {
+    // Get gather set cells
+    Range gather_set_cells;
+    rval = mb.get_entities_by_type(gather_set, MBPOLYGON, gather_set_cells);
+    CHECK_ERR(rval);
+    CHECK_EQUAL((size_t)642, gather_set_cells.size());
+    CHECK_EQUAL((size_t)1, gather_set_cells.psize());
+
+    // Check vorticity0 tag values on 4 gather set cells: first cell, two median cells, and last cell
+    EntityHandle cell_ents[] = {gather_set_cells[0], gather_set_cells[320],
+                                gather_set_cells[321], gather_set_cells[641]};
+    double vorticity0_val[4 * layers];
+    rval = mb.tag_get_data(vorticity_tag0, cell_ents, 4, vorticity0_val);
+    CHECK_ERR(rval);
+
+    // Only check first two layers
+    // Layer 0
+    CHECK_REAL_EQUAL(3.629994, vorticity0_val[0 * layers], eps);
+    CHECK_REAL_EQUAL(0.131688, vorticity0_val[1 * layers], eps);
+    CHECK_REAL_EQUAL(-0.554888, vorticity0_val[2 * layers], eps);
+    CHECK_REAL_EQUAL(-0.554888, vorticity0_val[3 * layers], eps);
+    // Layer 1
+    CHECK_REAL_EQUAL(3.629944, vorticity0_val[0 * layers + 1], eps);
+    CHECK_REAL_EQUAL(0.131686, vorticity0_val[1 * layers + 1], eps);
+    CHECK_REAL_EQUAL(-0.554881, vorticity0_val[2 * layers + 1], eps);
+    CHECK_REAL_EQUAL(-0.554881, vorticity0_val[3 * layers + 1], eps);
+  }
+#endif
 }
 
 void get_options(std::string& opts)

diff --git a/test/parallel/gcrm_par.cpp b/test/parallel/gcrm_par.cpp
index b2dae0a..c6b6fe4 100644
--- a/test/parallel/gcrm_par.cpp
+++ b/test/parallel/gcrm_par.cpp
@@ -54,8 +54,8 @@ int main(int argc, char* argv[])
   result += RUN_TEST(test_read_mesh_parallel_rcbzoltan);
 #endif
 
-  //result += RUN_TEST(test_gather_onevar_on_rank0);
-  //result += RUN_TEST(test_gather_onevar_on_rank1);
+  result += RUN_TEST(test_gather_onevar_on_rank0);
+  result += RUN_TEST(test_gather_onevar_on_rank1);
 
   result += RUN_TEST(test_multiple_loads_of_same_file);
 
@@ -508,33 +508,39 @@ void gather_one_cell_var(int gather_set_rank)
     assert(gather_set != 0);
   }
 
-  Tag ke_tag0, gid_tag;
-  rval = mb.tag_get_handle("ke0", 1, MB_TYPE_DOUBLE, ke_tag0, MB_TAG_DENSE);
+  Tag vorticity_tag0, gid_tag;
+  rval = mb.tag_get_handle("vorticity0", layers, MB_TYPE_DOUBLE, vorticity_tag0, MB_TAG_DENSE);
   CHECK_ERR(rval);
 
   rval = mb.tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gid_tag, MB_TAG_DENSE);
   CHECK_ERR(rval);
 
-  pcomm->gather_data(cells_owned, ke_tag0, gid_tag, gather_set, gather_set_rank);
+  pcomm->gather_data(cells_owned, vorticity_tag0, gid_tag, gather_set, gather_set_rank);
 
   if (gather_set_rank == rank) {
     // Get gather set cells
     Range gather_set_cells;
     rval = mb.get_entities_by_type(gather_set, MBPOLYGON, gather_set_cells);
     CHECK_EQUAL((size_t)642, gather_set_cells.size());
-    CHECK_EQUAL((size_t)2, gather_set_cells.psize());
-
-    // Check ke0 tag values on 4 gather set cells: first pentagon, last pentagon,
-    // first hexagon and last hexagon
-    EntityHandle cell_ents[] = {gather_set_cells[0], gather_set_cells[11],
-                                gather_set_cells[12], gather_set_cells[641]};
-    double ke0_val[4];
-    rval = mb.tag_get_data(ke_tag0, &cell_ents[0], 4, ke0_val);
-
-    CHECK_REAL_EQUAL(15.001, ke0_val[0], eps);
-    CHECK_REAL_EQUAL(15.012, ke0_val[1], eps);
-    CHECK_REAL_EQUAL(16.013, ke0_val[2], eps);
-    CHECK_REAL_EQUAL(16.642, ke0_val[3], eps);
+    CHECK_EQUAL((size_t)1, gather_set_cells.psize());
+
+    // Check vorticity0 tag values on 4 gather set cells: first cell, two median cells, and last cell
+    EntityHandle cell_ents[] = {gather_set_cells[0], gather_set_cells[320],
+                                gather_set_cells[321], gather_set_cells[641]};
+    double vorticity0_val[4 * layers];
+    rval = mb.tag_get_data(vorticity_tag0, &cell_ents[0], 4, vorticity0_val);
+
+    // Only check first two layers
+    // Layer 0
+    CHECK_REAL_EQUAL(3.629994, vorticity0_val[0 * layers], eps);
+    CHECK_REAL_EQUAL(0.131688, vorticity0_val[1 * layers], eps);
+    CHECK_REAL_EQUAL(-0.554888, vorticity0_val[2 * layers], eps);
+    CHECK_REAL_EQUAL(-0.554888, vorticity0_val[3 * layers], eps);
+    // Layer 1
+    CHECK_REAL_EQUAL(3.629944, vorticity0_val[0 * layers + 1], eps);
+    CHECK_REAL_EQUAL(0.131686, vorticity0_val[1 * layers + 1], eps);
+    CHECK_REAL_EQUAL(-0.554881, vorticity0_val[2 * layers + 1], eps);
+    CHECK_REAL_EQUAL(-0.554881, vorticity0_val[3 * layers + 1], eps);
   }
 }

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