[cig-commits] r21859 - in short/3D/PyLith/trunk: libsrc/pylith/meshio unittests/libtests/meshio

brad at geodynamics.org brad at geodynamics.org
Sat Apr 13 09:36:29 PDT 2013


Author: brad
Date: 2013-04-13 09:36:29 -0700 (Sat, 13 Apr 2013)
New Revision: 21859

Modified:
   short/3D/PyLith/trunk/libsrc/pylith/meshio/HDF5.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIO.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOAscii.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOCubit.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOLagrit.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFile.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFileAscii.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFileBinary.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/Xdmf.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.hh
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOAscii.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOCubit.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOLagrit.cc
Log:
Code cleanup. Update MeshIO to use PETSc DM.

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/HDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/HDF5.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/HDF5.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -20,6 +20,8 @@
 
 #include "HDF5.hh" // implementation of class methods
 
+#include "pylith/utils/petscerror.h" // USES PYLITH_METHOD_BEGIN/END
+
 #include <cstring> // USES strlen()
 #include <stdexcept> // USES std::runtime_error
 #include <sstream> // USES std::ostringstream
@@ -41,6 +43,8 @@
 pylith::meshio::HDF5::HDF5(const char* filename,
 			   hid_t mode)
 { // constructor
+  PYLITH_METHOD_BEGIN;
+
   if (H5F_ACC_TRUNC == mode) {
     _file = H5Fcreate(filename, mode, H5P_DEFAULT, H5P_DEFAULT);
     if (_file < 0) {
@@ -57,6 +61,8 @@
       throw std::runtime_error(msg.str());
     } // if
   } // if/else
+
+  PYLITH_METHOD_END;
 } // constructor
 
 // ----------------------------------------------------------------------
@@ -72,6 +78,8 @@
 pylith::meshio::HDF5::open(const char* filename,
 			   hid_t mode)
 { // open
+  PYLITH_METHOD_BEGIN;
+
   assert(filename);
 
   if (_file >= 0) {
@@ -94,6 +102,8 @@
       throw std::runtime_error(msg.str());
     } // if
   } // if/else
+
+  PYLITH_METHOD_END;
 } // constructor
 
 // ----------------------------------------------------------------------
@@ -101,12 +111,16 @@
 void
 pylith::meshio::HDF5::close(void)
 { // close
+  PYLITH_METHOD_BEGIN;
+
   if (_file >= 0) {
     herr_t err = H5Fclose(_file);
     if (err < 0) 
       throw std::runtime_error("Could not close HDF5 file.");
   } // if
   _file = -1;
+
+  PYLITH_METHOD_END;
 } // close
 
 // ----------------------------------------------------------------------
@@ -122,6 +136,8 @@
 bool
 pylith::meshio::HDF5::hasGroup(const char* name)
 { // hasGroup
+  PYLITH_METHOD_BEGIN;
+
   assert(isOpen());
   assert(name);
 
@@ -138,7 +154,7 @@
     assert(err >= 0);
   } // if
   
-  return exists;
+  PYLITH_METHOD_RETURN(exists);
 } // hasGroup
 
 // ----------------------------------------------------------------------
@@ -146,6 +162,8 @@
 bool
 pylith::meshio::HDF5::hasDataset(const char* name)
 { // hasDataset
+  PYLITH_METHOD_BEGIN;
+
   assert(isOpen());
   assert(name);
 
@@ -162,7 +180,7 @@
     assert(err >= 0);
   } // if
   
-  return exists;
+  PYLITH_METHOD_RETURN(exists);
 } // hasDataset
 
 // ----------------------------------------------------------------------
@@ -173,6 +191,8 @@
 				     const char* parent,
 				     const char* name)
 { // getDatasetDims
+  PYLITH_METHOD_BEGIN;
+
   assert(dims);
   assert(ndims);
   assert(isOpen());
@@ -228,6 +248,8 @@
 	<< parent << "/" << name << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch  
+
+  PYLITH_METHOD_END;
 } // getDatasetDims
 
 // ----------------------------------------------------------------------
@@ -236,6 +258,8 @@
 pylith::meshio::HDF5::getGroupDatasets(string_vector* names,
 				       const char* parent)
 { // getGroupDatasets
+  PYLITH_METHOD_BEGIN;
+
   assert(names);
   assert(isOpen());
 
@@ -280,6 +304,8 @@
 	<< parent << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch  
+
+  PYLITH_METHOD_END;
 } // getGroupDatasets
 
 // ----------------------------------------------------------------------
@@ -287,6 +313,8 @@
 void
 pylith::meshio::HDF5::createGroup(const char* name)
 { // createGroup
+  PYLITH_METHOD_BEGIN;
+
   assert(name);
 
 #if defined(PYLITH_HDF5_USE_API_18)
@@ -306,6 +334,8 @@
     msg << "Could not close group '" << name << "'.";
     throw std::runtime_error(msg.str());
   } // if
+
+  PYLITH_METHOD_END;
 } // createGroup
 
 // ----------------------------------------------------------------------
@@ -316,6 +346,8 @@
 				     const void* value,
 				     hid_t datatype)
 { // writeAttribute
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(value);
@@ -364,6 +396,8 @@
     msg << err.what() << " attribute '" << name << "' of '" << parent << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // writeAttribute
 
 // ----------------------------------------------------------------------
@@ -375,6 +409,8 @@
 				     const void* value,
 				     hid_t datatype)
 { // writeAttribute
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(value);
@@ -423,6 +459,8 @@
     msg << err.what() << " attribute '" << name << "' of '" << parent << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // writeAttribute
 
 // ----------------------------------------------------------------------
@@ -433,6 +471,8 @@
 				    void* value,
 				    hid_t datatype)
 { // readAttribute
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(value);
@@ -475,6 +515,8 @@
     msg << err.what() << " attribute '" << name << "' of '" << parent << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // readAttribute
 
 // ----------------------------------------------------------------------
@@ -484,6 +526,8 @@
 				     const char* name,
 				     const char* value)
 { // writeAttribute
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(value);
@@ -544,6 +588,8 @@
     msg << err.what() << " attribute '" << name << "' of '" << parent << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // writeAttribute
 
 // ----------------------------------------------------------------------
@@ -554,6 +600,8 @@
 				     const char* name,
 				     const char* value)
 { // writeAttribute
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(value);
@@ -614,6 +662,8 @@
     msg << err.what() << " attribute '" << name << "' of '" << parent << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // writeAttribute
 
 // ----------------------------------------------------------------------
@@ -622,6 +672,8 @@
 pylith::meshio::HDF5::readAttribute(const char* parent,
 				    const char* name)
 { // readAttribute
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
 
@@ -677,7 +729,7 @@
     throw std::runtime_error(msg.str());
   } // try/catch
 
-  return std::string(value);
+  PYLITH_METHOD_RETURN(std::string(value));
 } // readAttribute
 
 // ----------------------------------------------------------------------
@@ -690,6 +742,8 @@
 				    const int ndims,
 				    hid_t datatype)
 { // createDataset
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(maxDims);
@@ -766,6 +820,8 @@
     msg << "Unknown  occurred while creating dataset '" << name << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // createDataset
 
 // ----------------------------------------------------------------------
@@ -780,6 +836,8 @@
 					const int chunk,
 					hid_t datatype)
 { // writeDatasetSlice
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(data);
@@ -874,6 +932,8 @@
 	<< parent << "/" << name << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // writeDatasetSlice
 
 // ----------------------------------------------------------------------
@@ -887,6 +947,8 @@
 				       const int chunk,
 				       hid_t datatype)
 { // readDatasetSlice
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(data);
@@ -987,6 +1049,8 @@
 	<< parent << "/" << name << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // readDatasetSlice
 
 // ----------------------------------------------------------------------
@@ -1000,6 +1064,8 @@
 					       const int ndims,
 					       hid_t datatype)
 { // createDatasetRawExternal
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(filename);
@@ -1075,6 +1141,8 @@
     msg << "Unknown  occurred while creating dataset '" << name << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // createDatasetRawExternal
 
 // ----------------------------------------------------------------------
@@ -1086,6 +1154,8 @@
 					       const hsize_t* dims,
 					       const int ndims)
 { // extendDatasetRawExternal
+  PYLITH_METHOD_BEGIN;
+
   assert(parent);
   assert(name);
   assert(dims);
@@ -1136,6 +1206,8 @@
     msg << "Unknown  occurred while updating dataset '" << name << "'.";
     throw std::runtime_error(msg.str());
   } // try/catch
+
+  PYLITH_METHOD_END;
 } // extendDatasetRawExternal
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -53,11 +53,11 @@
 
   assert(mesh);
   assert(coordinates);
-  MPI_Comm       comm     = mesh->comm();
-  PetscInt       dim      = meshDim;
-  PetscMPIInt    commRank = mesh->commRank();
-  const bool     useSieve = true;
+  MPI_Comm comm  = mesh->comm();
+  PetscInt dim  = meshDim;
+  PetscMPIInt commRank = mesh->commRank();
   PetscErrorCode err;
+  const bool useSieve = true; // :TODO: Remove sieve
 
   { // Check to make sure every vertex is in at least one cell.
     // This is required by Sieve
@@ -79,94 +79,82 @@
 
   /* Sieve */
   if (useSieve) {
-  MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
-  mesh->createSieveMesh(dim);
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
-  assert(!sieveMesh.isNull());
+    MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
+    mesh->createSieveMesh(dim);
+    const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+    assert(!sieveMesh.isNull());
 
-  ALE::Obj<SieveMesh::sieve_type> sieve = 
-    new SieveMesh::sieve_type(mesh->comm());
-  sieveMesh->setSieve(sieve);
+    ALE::Obj<SieveMesh::sieve_type> sieve = new SieveMesh::sieve_type(mesh->comm());
+    sieveMesh->setSieve(sieve);
 
-  if (0 == commRank || isParallel) {
-    assert(coordinates->size() == numVertices*spaceDim);
-    assert(cells.size() == numCells*numCorners);
-    if (!interpolate) {
-      // Create the ISieve
-      sieve->setChart(SieveMesh::sieve_type::chart_type(0, 
-							numCells+numVertices));
-      // Set cone and support sizes
-      for (int c = 0; c < numCells; ++c)
-	sieve->setConeSize(c, numCorners);
-      sieve->symmetrizeSizes(numCells, numCorners, 
-			     const_cast<int*>(&cells[0]), numCells);
-      // Allocate point storage
-      sieve->allocate();
-      // Fill up cones
-      int *cone  = new int[numCorners];
-      int *coneO = new int[numCorners];
-      for (int v = 0; v < numCorners; ++v)
-	coneO[v] = 1;
-      for (int c = 0; c < numCells; ++c) {
-        for (int v = 0; v < numCorners; ++v)
-	  cone[v] = cells[c*numCorners+v]+numCells;
-        sieve->setCone(cone, c);
-        sieve->setConeOrientation(coneO, c);
-      } // for
-      delete[] cone; cone = 0;
-      delete[] coneO; coneO = 0;
-      // Symmetrize to fill up supports
-      sieve->symmetrize();
-    } else {
-      // Same old thing
-      ALE::Obj<SieveFlexMesh::sieve_type> s =
-	new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
-      ALE::Obj<SieveFlexMesh::arrow_section_type> orientation = new SieveFlexMesh::arrow_section_type(sieve->comm(), sieve->debug());
+    if (0 == commRank || isParallel) {
+      assert(coordinates->size() == numVertices*spaceDim);
+      assert(cells.size() == numCells*numCorners);
+      if (!interpolate) {
+	// Create the ISieve
+	sieve->setChart(SieveMesh::sieve_type::chart_type(0, numCells+numVertices));
+	// Set cone and support sizes
+	for (int c = 0; c < numCells; ++c)
+	  sieve->setConeSize(c, numCorners);
+	sieve->symmetrizeSizes(numCells, numCorners, const_cast<int*>(&cells[0]), numCells);
+	// Allocate point storage
+	sieve->allocate();
+	// Fill up cones
+	int *cone  = new int[numCorners];
+	int *coneO = new int[numCorners];
+	for (int v = 0; v < numCorners; ++v)
+	  coneO[v] = 1;
+	for (int c = 0; c < numCells; ++c) {
+	  for (int v = 0; v < numCorners; ++v)
+	    cone[v] = cells[c*numCorners+v]+numCells;
+	  sieve->setCone(cone, c);
+	  sieve->setConeOrientation(coneO, c);
+	} // for
+	delete[] cone; cone = 0;
+	delete[] coneO; coneO = 0;
+	// Symmetrize to fill up supports
+	sieve->symmetrize();
+      } else {
+	// Same old thing
+	ALE::Obj<SieveFlexMesh::sieve_type> s = new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
+	ALE::Obj<SieveFlexMesh::arrow_section_type> orientation = new SieveFlexMesh::arrow_section_type(sieve->comm(), sieve->debug());
 
-      s->setDebug(2);
-      ALE::SieveBuilder<SieveFlexMesh>::buildTopology(s, meshDim, 
-                                                      numCells, 
-                                                      const_cast<int*>(&cells[0]), 
-                                                      numVertices, 
-                                                      interpolate,
-                                                      numCorners,
-                                                      -1,
-                                                      orientation);
-      std::map<SieveFlexMesh::point_type,SieveFlexMesh::point_type> renumbering;
-      ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
-    } // if/else
-    if (!interpolate) {
-      // Optimized stratification
-      const ALE::Obj<SieveMesh::label_type>& height = 
-	sieveMesh->createLabel("height");
-      const ALE::Obj<SieveMesh::label_type>& depth =
-	sieveMesh->createLabel("depth");
-
-      for(int c = 0; c < numCells; ++c) {
-        height->setCone(0, c);
-        depth->setCone(1, c);
-      } // for
-      for(int v = numCells; v < numCells+numVertices; ++v) {
-        height->setCone(1, v);
-        depth->setCone(0, v);
-      } // for
-      sieveMesh->setHeight(1);
-      sieveMesh->setDepth(1);
+	s->setDebug(2);
+	ALE::SieveBuilder<SieveFlexMesh>::buildTopology(s, meshDim, numCells, const_cast<int*>(&cells[0]), numVertices, interpolate,
+							numCorners, -1, orientation);
+	std::map<SieveFlexMesh::point_type,SieveFlexMesh::point_type> renumbering;
+	ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
+      } // if/else
+      if (!interpolate) {
+	// Optimized stratification
+	const ALE::Obj<SieveMesh::label_type>& height = sieveMesh->createLabel("height");
+	const ALE::Obj<SieveMesh::label_type>& depth = sieveMesh->createLabel("depth");
+	
+	for(int c = 0; c < numCells; ++c) {
+	  height->setCone(0, c);
+	  depth->setCone(1, c);
+	} // for
+	for(int v = numCells; v < numCells+numVertices; ++v) {
+	  height->setCone(1, v);
+	  depth->setCone(0, v);
+	} // for
+	sieveMesh->setHeight(1);
+	sieveMesh->setDepth(1);
+      } else {
+	sieveMesh->stratify();
+      } // if/else
     } else {
+      sieveMesh->getSieve()->setChart(SieveMesh::sieve_type::chart_type());
+      sieveMesh->getSieve()->allocate();
       sieveMesh->stratify();
     } // if/else
-  } else {
-    sieveMesh->getSieve()->setChart(SieveMesh::sieve_type::chart_type());
-    sieveMesh->getSieve()->allocate();
-    sieveMesh->stratify();
-  } // if/else
 
-  ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim, &(*coordinates)[0]);
-  sieveMesh->getFactory()->clear();
-  }
+    ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim, &(*coordinates)[0]);
+    sieveMesh->getFactory()->clear();
+  } // if
 
   /* DMPlex */
-  DM        dmMesh;
+  PetscDM dmMesh;
   PetscBool pInterpolate = interpolate ? PETSC_TRUE : PETSC_FALSE;
 
   err = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, pInterpolate, &cells[0], spaceDim, &(*coordinates)[0], &dmMesh);CHECK_PETSC_ERROR(err);

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIO.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIO.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -21,16 +21,19 @@
 #include "MeshIO.hh" // implementation of class methods
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/Stratum.hh" // USES Stratum
+
 #include "pylith/utils/array.hh" // USES scalar_array, int_array
+#include "pylith/utils/petscerror.h" // USES PYLITH_METHOD_BEGIN/END
 
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+
 #include "Selection.hh" // USES boundary()
 
 #include <cassert> // USES assert()
 #include <sstream> // USES std::ostringstream
 #include <stdexcept> // USES std::runtime_error
 
-#include <utils/petscerror.h> // USES CHECK_PETSC_ERROR
-
 // ----------------------------------------------------------------------
 typedef pylith::topology::Mesh::SieveMesh SieveMesh;
 typedef pylith::topology::Mesh::RealSection RealSection;
@@ -64,7 +67,7 @@
 int
 pylith::meshio::MeshIO::getMeshDim(void) const
 { // getMeshDim
-  return (0 != _mesh) ? _mesh->dimension() : 0;
+  return (_mesh) ? _mesh->dimension() : 0;
 } // getMeshDim
 
 // ----------------------------------------------------------------------
@@ -74,6 +77,7 @@
 { // read
   PYLITH_METHOD_BEGIN;
 
+  assert(mesh);
   assert(!_mesh);
 
   _mesh = mesh;
@@ -90,11 +94,16 @@
 void 
 pylith::meshio::MeshIO::write(topology::Mesh* const mesh)
 { // write
-  assert(0 == _mesh);
+  PYLITH_METHOD_BEGIN;
 
+  assert(mesh);
+  assert(!_mesh);
+
   _mesh = mesh;
   _write();
   _mesh = 0;
+
+  PYLITH_METHOD_END;
 } // write
 
 // ----------------------------------------------------------------------
@@ -104,21 +113,21 @@
 				     int* numVertices,
 				     int* spaceDim) const
 { // _getVertices
-  assert(0 != coordinates);
-  assert(0 != numVertices);
-  assert(0 != spaceDim);
-  assert(0 != _mesh);
+  PYLITH_METHOD_BEGIN;
 
+  assert(coordinates);
+  assert(numVertices);
+  assert(spaceDim);
+  assert(_mesh);
+
+#if 0 // Sieve
   const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
   assert(!sieveMesh.isNull());
 
-  const ALE::Obj<SieveMesh::label_sequence>& vertices = 
-    sieveMesh->depthStratum(0);
+  const ALE::Obj<SieveMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
   assert(!vertices.isNull());
-  const SieveMesh::label_sequence::iterator verticesBegin = 
-    vertices->begin();
-  const SieveMesh::label_sequence::iterator verticesEnd = 
-    vertices->end();
+  const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+  const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
   const ALE::Obj<RealSection>& coordsField =
     sieveMesh->hasRealSection("coordinates_dimensioned") ?
     sieveMesh->getRealSection("coordinates_dimensioned") :
@@ -143,6 +152,35 @@
     for (int iDim=0; iDim < *spaceDim; ++iDim)
       (*coordinates)[i++] = vertexCoords[iDim];
   } // for
+#endif
+
+  const spatialdata::geocoords::CoordSys* cs = _mesh->coordsys();assert(cs);
+  *spaceDim = cs->spaceDim();
+
+  PetscDM dmMesh = _mesh->dmMesh();assert(dmMesh);
+  PetscVec coordVec = NULL;
+  PetscScalar* coordArray = NULL;
+  PetscInt coordSize = 0;    
+  PylithScalar lengthScale = 1.0;
+  PetscErrorCode err = 0;
+
+  // Get length scale for dimensioning
+  err = DMPlexGetScale(dmMesh, PETSC_UNIT_LENGTH, &lengthScale);
+
+  // Get coordinates and dimensionalize values
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
+  err = VecGetLocalSize(coordVec, &coordSize);CHECK_PETSC_ERROR(err);
+  assert(coordSize % *spaceDim == 0);
+  *numVertices = coordSize / *spaceDim;
+
+  coordinates->resize(coordSize);
+  for (PetscInt i=0; i < coordSize; ++i) {
+    (*coordinates)[i] = coordArray[i]*lengthScale;
+  } // for
+  err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
+
+  PYLITH_METHOD_END;
 } // _getVertices
 
 // ----------------------------------------------------------------------
@@ -153,24 +191,23 @@
 				  int* numCorners,
 				  int* meshDim) const
 { // _getCells
+  PYLITH_METHOD_BEGIN;
 
-  assert(0 != cells);
-  assert(0 != numCells);
-  assert(0 != meshDim);
-  assert(0 != _mesh);
+  assert(cells);
+  assert(numCells);
+  assert(meshDim);
+  assert(_mesh);
 
+#if 0 // Sieve
   const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
   assert(!sieveMesh.isNull());
 
   const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
   assert(!sieve.isNull());
-  const ALE::Obj<SieveMesh::label_sequence>& meshCells = 
-    sieveMesh->heightStratum(0);
+  const ALE::Obj<SieveMesh::label_sequence>& meshCells = sieveMesh->heightStratum(0);
   assert(!meshCells.isNull());
-  const SieveMesh::label_sequence::iterator meshCellsBegin = 
-    meshCells->begin();
-  const SieveMesh::label_sequence::iterator meshCellsEnd = 
-    meshCells->end();
+  const SieveMesh::label_sequence::iterator meshCellsBegin = meshCells->begin();
+  const SieveMesh::label_sequence::iterator meshCellsEnd = meshCells->end();
 
   *meshDim = _mesh->dimension();
   *numCells = meshCells->size();
@@ -197,6 +234,44 @@
     }
     pV.clear();
   } // for
+#endif // Sieve
+
+  PetscDM dmMesh = _mesh->dmMesh();assert(dmMesh);
+  topology::Stratum cellsStratum(dmMesh, topology::Stratum::HEIGHT, 0);
+  const PetscInt cStart = cellsStratum.begin();
+  const PetscInt cEnd = cellsStratum.end();
+  
+  topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
+  
+  *numCells = _mesh->numCells();
+  *numCorners = _mesh->coneSize();
+  *meshDim = _mesh->dimension();
+  assert(cellsStratum.size() == *numCells);
+
+  cells->resize((*numCells)*(*numCorners));
+
+  PetscIS globalVertexNumbers = NULL;
+  const PetscInt* gvertex = NULL;
+  const PetscInt* cone = NULL;
+  PetscInt coneSize = 0, v = 0, count = 0;
+  PetscErrorCode err = 0;
+
+  err = DMPlexGetVertexNumbering(dmMesh, &globalVertexNumbers);CHECK_PETSC_ERROR(err);
+  err = ISGetIndices(globalVertexNumbers, &gvertex);CHECK_PETSC_ERROR(err);
+  for (PetscInt c = cStart, index = 0; c < cEnd; ++c) {
+    err = DMPlexGetConeSize(dmMesh, c, &coneSize);CHECK_PETSC_ERROR(err);
+    assert(coneSize == *numCorners);
+
+    err = DMPlexGetCone(dmMesh, c, &cone);CHECK_PETSC_ERROR(err);
+    for(PetscInt p = 0; p < coneSize; ++p) {
+      const PetscInt gv = gvertex[cone[p]-vStart];
+      (*cells)[index++] = gv < 0 ? -(gv+1) : gv;
+    } // for
+  } // for  
+
+  PYLITH_METHOD_END;
 } // _getCells
 
 // ----------------------------------------------------------------------
@@ -204,13 +279,13 @@
 void
 pylith::meshio::MeshIO::_setMaterials(const int_array& materialIds)
 { // _setMaterials
-  assert(0 != _mesh);
+  PYLITH_METHOD_BEGIN;
+
+  assert(_mesh);
   PetscErrorCode err;
 
-  ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
-  ///logger.setDebug(2);
-  logger.stagePush("MeshLabels");
   if (!_mesh->commRank()) {
+#if 1 // Sieve
     const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
     if (!sieveMesh.isNull()) {
       const ALE::Obj<SieveMesh::label_type>& labelMaterials = sieveMesh->createLabel("material-id");
@@ -231,24 +306,26 @@
         sieveMesh->setValue(labelMaterials, *e_iter, materialIds[i++]);
       }
     }
+#endif // Sieve
 
-    DM       dmMesh = _mesh->dmMesh();
-    PetscInt cStart, cEnd;
+    PetscDM dmMesh = _mesh->dmMesh();assert(dmMesh);
+    topology::Stratum cellsStratum(dmMesh, topology::Stratum::HEIGHT, 0);
+    const PetscInt cStart = cellsStratum.begin();
+    const PetscInt cEnd = cellsStratum.end();
 
-    assert(dmMesh);
-    err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
-    if ((cEnd - cStart) != materialIds.size()) {
+    if (cellsStratum.size() != materialIds.size()) {
       std::ostringstream msg;
       msg << "Mismatch in size of materials identifier array ("
           << materialIds.size() << ") and number of cells in mesh ("<< (cEnd - cStart) << ").";
       throw std::runtime_error(msg.str());
     } // if
+    PetscErrorCode err = 0;
     for(PetscInt c = cStart; c < cEnd; ++c) {
       err = DMPlexSetLabelValue(dmMesh, "material-id", c, materialIds[c-cStart]);CHECK_PETSC_ERROR(err);
-    }
+    } // for
   } // if
-  logger.stagePop();
-  ///logger.setDebug(0);
+
+  PYLITH_METHOD_END;
 } // _setMaterials
 
 // ----------------------------------------------------------------------
@@ -256,9 +333,12 @@
 void
 pylith::meshio::MeshIO::_getMaterials(int_array* materialIds) const
 { // _getMaterials
-  assert(0 != materialIds);
-  assert(0 != _mesh);
+  PYLITH_METHOD_BEGIN;
 
+  assert(materialIds);
+  assert(_mesh);
+
+#if 0 // Sieve
   const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
   assert(!sieveMesh.isNull());
 
@@ -282,6 +362,22 @@
       ++e_iter)
     (*materialIds)[i++] = 
       sieveMesh->getValue(labelMaterials, *e_iter, idDefault);
+#endif // Sieve
+
+  PetscDM dmMesh = _mesh->dmMesh();assert(dmMesh);
+  topology::Stratum cellsStratum(dmMesh, topology::Stratum::HEIGHT, 0);
+  const PetscInt cStart = cellsStratum.begin();
+  const PetscInt cEnd = cellsStratum.end();
+
+  materialIds->resize(cellsStratum.size());
+  PetscErrorCode err = 0;
+  PetscInt matId = 0;
+  for(PetscInt c = cStart, index = 0; c < cEnd; ++c) {
+    err = DMPlexGetLabelValue(dmMesh, "material-id", c, &matId);CHECK_PETSC_ERROR(err);
+    (*materialIds)[index++] = matId;
+  } // for
+
+  PYLITH_METHOD_END;
 } // _getMaterials
 
 // ----------------------------------------------------------------------
@@ -291,11 +387,11 @@
 				  const GroupPtType type,
 				  const int_array& points)
 { // _setGroup
-  assert(0 != _mesh);
-  PetscErrorCode err;
-  ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
-  logger.stagePush("MeshIntSections");
+  PYLITH_METHOD_BEGIN;
 
+  assert(_mesh);
+
+#if 1 // Sieve
   const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
   if (!sieveMesh.isNull()) {
     if (sieveMesh->hasIntSection(name)) {
@@ -322,15 +418,16 @@
     } // if/else
     sieveMesh->allocate(groupField);
   }
+#endif // Sieve
 
-  DM             dmMesh    = _mesh->dmMesh();
+  PetscDM dmMesh = _mesh->dmMesh();assert(dmMesh);
   const PetscInt numPoints = points.size();
+  PetscErrorCode err;
 
-  assert(dmMesh);
   if (CELL == type) {
     for(PetscInt p = 0; p < numPoints; ++p) {
       err = DMPlexSetLabelValue(dmMesh, name.c_str(), points[p], 1);CHECK_PETSC_ERROR(err);
-    }
+    } // for
   } else if (VERTEX == type) {
     PetscInt cStart, cEnd, numCells;
 
@@ -338,9 +435,10 @@
     numCells = cEnd - cStart;
     for(PetscInt p = 0; p < numPoints; ++p) {
       err = DMPlexSetLabelValue(dmMesh, name.c_str(), numCells+points[p], 1);CHECK_PETSC_ERROR(err);
-    }
-  }
-  logger.stagePop();
+    } // for
+  } // if/else
+
+  PYLITH_METHOD_END;
 } // _setGroup
 
 // ----------------------------------------------------------------------
@@ -348,9 +446,13 @@
 void
 pylith::meshio::MeshIO::_distributeGroups()
 { // _distributeGroups
-  assert(0 != _mesh);
-  /* dmMesh does not need to broadcast the label names. They come from proc 0 */
+  PYLITH_METHOD_BEGIN;
 
+  assert(_mesh);
+
+  // dmMesh does not need to broadcast the label names. They come from proc 0.
+
+#if 1 // Sieve
   const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
   if (!sieveMesh.isNull()) {
     if (!sieveMesh->commRank()) {
@@ -388,6 +490,9 @@
       } // for
     } // if/else
   }
+#endif // Sieve
+
+  PYLITH_METHOD_END;
 } // _distributeGroups
 
 // ----------------------------------------------------------------------
@@ -395,9 +500,12 @@
 void
 pylith::meshio::MeshIO::_getGroupNames(string_vector* names) const
 { // _getGroups
-  assert(0 != names);
-  assert(0 != _mesh);
+  PYLITH_METHOD_BEGIN;
 
+  assert(names);
+  assert(_mesh);
+
+#if 0 // Sieve
   const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
   assert(!sieveMesh.isNull());
 
@@ -412,19 +520,39 @@
        name != namesEnd;
        ++name)
     (*names)[i++] = *name;
+#endif // Sieve
+
+  PetscDM dmMesh = _mesh->dmMesh();assert(dmMesh);
+  PetscInt numGroups = 0;
+  PetscErrorCode err = 0;
+  err = DMPlexGetNumLabels(dmMesh, &numGroups);CHECK_PETSC_ERROR(err);
+  numGroups -= 2; // Remove depth and material labels.
+  names->resize(numGroups);
+
+  for (int iGroup=0, iLabel=numGroups-1; iGroup < numGroups; ++iGroup, --iLabel) {
+    const char* namestr = NULL;
+    err = DMPlexGetLabelName(dmMesh, iLabel, &namestr);CHECK_PETSC_ERROR(err);
+    (*names)[iGroup] = namestr;
+    std::cout << "GROUP " << iGroup << ": " << (*names)[iGroup] << std::endl;
+  } // for
+
+  PYLITH_METHOD_END;
 } // _getGroups
 
 // ----------------------------------------------------------------------
 // Get group entities
 void
 pylith::meshio::MeshIO::_getGroup(int_array* points,
-				  GroupPtType* type,
+				  GroupPtType* groupType,
 				  const char *name) const
 { // _getGroup
-  assert(0 != points);
-  assert(0 != type);
-  assert(0 != _mesh);
+  PYLITH_METHOD_BEGIN;
 
+  assert(points);
+  assert(groupType);
+  assert(_mesh);
+
+#if 0
   const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
   assert(!sieveMesh.isNull());
 
@@ -451,11 +579,11 @@
   ALE::Obj<SieveMesh::numbering_type> numbering;
 
   if (sieveMesh->height(firstPoint) == 0) {
-    *type = CELL;
+    *groupType = CELL;
     numbering = sieveMesh->getFactory()->getNumbering(sieveMesh, 
 						      sieveMesh->depth());
   } else {
-    *type = VERTEX;
+    *groupType = VERTEX;
     numbering = sieveMesh->getFactory()->getNumbering(sieveMesh, 0);
   } // if/else
   const int numPoints = groupField->size();
@@ -469,6 +597,43 @@
     if (groupField->getFiberDimension(*c_iter)) (*points)[i++] = 
       numbering->getIndex(*c_iter);
   } // for
+#endif
+
+  PetscDM dmMesh = _mesh->dmMesh();assert(dmMesh);
+  topology::Stratum cellsStratum(dmMesh, topology::Stratum::HEIGHT, 0);
+  const PetscInt cStart = cellsStratum.begin();
+  const PetscInt cEnd = cellsStratum.end();
+  const PetscInt numCells = cellsStratum.size();
+
+  PetscInt pStart, pEnd, firstPoint = 0;
+  PetscErrorCode err = 0;
+  err = DMPlexGetChart(dmMesh, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+  for(PetscInt p = pStart; p < pEnd; ++p) {
+    PetscInt val;
+    err = DMPlexGetLabelValue(dmMesh, name, p, &val);CHECK_PETSC_ERROR(err);
+    if (val >= 0) {
+      firstPoint = p;
+      break;
+    } // if
+  } // for
+  *groupType = (firstPoint >= cStart && firstPoint < cEnd) ? CELL : VERTEX;
+
+  PetscInt groupSize;
+  err = DMPlexGetStratumSize(dmMesh, name, 1, &groupSize);CHECK_PETSC_ERROR(err);
+  points->resize(groupSize);
+
+  const PetscInt offset = (VERTEX == *groupType) ? numCells : 0;
+  PetscIS groupIS = NULL;
+  const PetscInt* groupIndices = NULL;
+  err = DMPlexGetStratumIS(dmMesh, name, 1, &groupIS);CHECK_PETSC_ERROR(err);
+  err = ISGetIndices(groupIS, &groupIndices);CHECK_PETSC_ERROR(err);
+  for(PetscInt p = 0; p < groupSize; ++p) {
+    (*points)[p] = groupIndices[p]-offset;
+  } // for
+  err = ISRestoreIndices(groupIS, &groupIndices);CHECK_PETSC_ERROR(err);
+  err = ISDestroy(&groupIS);CHECK_PETSC_ERROR(err);
+
+  PYLITH_METHOD_END;
 } // _getGroup
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOAscii.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOAscii.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOAscii.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -36,7 +36,7 @@
 #include <sstream> // USES std::ostringstream
 
 // ----------------------------------------------------------------------
-const char* pylith::meshio::MeshIOAscii::groupTypeNames[] = {
+const char* pylith::meshio::MeshIOAscii::groupTypeNames[2] = {
   "vertices",
   "cells",
 };
@@ -61,7 +61,11 @@
 void
 pylith::meshio::MeshIOAscii::deallocate(void)
 { // deallocate
+  PYLITH_METHOD_BEGIN;
+
   MeshIO::deallocate();
+
+  PYLITH_METHOD_END;
 } // deallocate
   
 // ----------------------------------------------------------------------
@@ -69,6 +73,8 @@
 void
 pylith::meshio::MeshIOAscii::_read(void)
 { // _read
+  PYLITH_METHOD_BEGIN;
+
   const int commRank = _mesh->commRank();
   int meshDim = 0;
   int spaceDim = 0;
@@ -97,7 +103,7 @@
     
     buffer.str(parser.next());
     buffer >> token;
-    if (0 != strcasecmp(token.c_str(), "mesh")) {
+    if (strcasecmp(token.c_str(), "mesh")) {
       std::ostringstream msg;
       msg << "Expected 'mesh' token but encountered '" << token << "'\n";
       throw std::runtime_error(msg.str());
@@ -162,8 +168,7 @@
 	buffer >> token;
       } // while
       if (token != "}")
-	throw std::runtime_error("I/O error occurred while parsing mesh " \
-				 "tokens.");
+	throw std::runtime_error("I/O error occurred while parsing mesh tokens.");
     } catch (const std::exception& err) {
       std::ostringstream msg;
       msg << "Error occurred while reading PyLith mesh ASCII file '"
@@ -178,12 +183,12 @@
     } // catch
     filein.close();
   } else {
-    MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
-			   cells, numCells, numCorners, meshDim,
-			   _interpolate);
+    MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim, cells, numCells, numCorners, meshDim, _interpolate);
     _setMaterials(materialIds);
   } // if/else
   _distributeGroups();
+
+  PYLITH_METHOD_END;
 } // read
 
 // ----------------------------------------------------------------------
@@ -191,6 +196,8 @@
 void
 pylith::meshio::MeshIOAscii::_write(void) const
 { // write
+  PYLITH_METHOD_BEGIN;
+
   std::ofstream fileout(_filename.c_str());
   if (!fileout.is_open() || !fileout.good()) {
     std::ostringstream msg;
@@ -215,6 +222,8 @@
 
   fileout << "}\n";
   fileout.close();
+
+  PYLITH_METHOD_END;
 } // write
 
 // ----------------------------------------------------------------------
@@ -225,10 +234,12 @@
 					   int* numVertices, 
 					   int* numDims) const
 { // _readVertices
-  assert(0 != coordinates);
-  assert(0 != numVertices);
-  assert(0 != numDims);
+  PYLITH_METHOD_BEGIN;
 
+  assert(coordinates);
+  assert(numVertices);
+  assert(numDims);
+
   std::string token;
   std::istringstream buffer;
   const int maxIgnore = 1024;
@@ -270,6 +281,8 @@
   } // while
   if (token != "}")
     throw std::runtime_error("I/O error while parsing vertices.");
+
+  PYLITH_METHOD_END;
 } // _readVertices
 
 // ----------------------------------------------------------------------
@@ -277,6 +290,8 @@
 void
 pylith::meshio::MeshIOAscii::_writeVertices(std::ostream& fileout) const
 { // _writeVertices
+  PYLITH_METHOD_BEGIN;
+
   int spaceDim = 0;
   int numVertices = 0;
   scalar_array coordinates;
@@ -300,6 +315,8 @@
   fileout
     << "    }\n"
     << "  }\n";
+
+  PYLITH_METHOD_END;
 } // _writeVertices
   
 // ----------------------------------------------------------------------
@@ -311,11 +328,13 @@
 					int* numCells, 
 					int* numCorners) const
 { // _readCells
-  assert(0 != cells);
-  assert(0 != materialIds);
-  assert(0 != numCells);
-  assert(0 != numCorners);
+  PYLITH_METHOD_BEGIN;
 
+  assert(cells);
+  assert(materialIds);
+  assert(numCells);
+  assert(numCorners);
+
   int dimension = 0;
 
   std::string token;
@@ -388,6 +407,8 @@
     materialIds->resize(size);
     (*materialIds) = 0;
   } // if
+
+  PYLITH_METHOD_END;
 } // _readCells
 
 // ----------------------------------------------------------------------
@@ -395,6 +416,8 @@
 void
 pylith::meshio::MeshIOAscii::_writeCells(std::ostream& fileout) const
 { // _writeCells
+  PYLITH_METHOD_BEGIN;
+
   int meshDim = 0;
   int numCells = 0;
   int numCorners = 0;
@@ -427,6 +450,8 @@
   fileout << "    }\n";  
 
   fileout << "  }\n";
+
+  PYLITH_METHOD_END;
 } // _writeCells
 
 // ----------------------------------------------------------------------
@@ -437,10 +462,12 @@
 					GroupPtType* type,
 					std::string* name) const
 { // _readGroup
-  assert(0 != points);
-  assert(0 != type);
-  assert(0 != name);
+  PYLITH_METHOD_BEGIN;
 
+  assert(points);
+  assert(type);
+  assert(name);
+
   std::string token;
   std::istringstream buffer;
   const int maxIgnore = 1024;
@@ -507,6 +534,8 @@
 
   if (!_useIndexZero)
     *points -= 1;
+
+  PYLITH_METHOD_END;
 } // _readGroup
 
 // ----------------------------------------------------------------------
@@ -515,6 +544,8 @@
 pylith::meshio::MeshIOAscii::_writeGroup(std::ostream& fileout,
 					 const char* name) const
 { // _writeGroup
+  PYLITH_METHOD_BEGIN;
+
   int_array points;
   GroupPtType type;
   _getGroup(&points, &type, name);
@@ -534,6 +565,8 @@
   fileout
     << "    }\n"
     << "  }\n";
+
+  PYLITH_METHOD_END;
 } // _writeGroup
   
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOCubit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOCubit.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOCubit.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -52,7 +52,11 @@
 void
 pylith::meshio::MeshIOCubit::deallocate(void)
 { // deallocate
+  PYLITH_METHOD_BEGIN;
+
   MeshIO::deallocate();
+
+  PYLITH_METHOD_END;
 } // deallocate
   
 // ----------------------------------------------------------------------
@@ -60,6 +64,8 @@
 void
 pylith::meshio::MeshIOCubit::_read(void)
 { // _read
+  PYLITH_METHOD_BEGIN;
+
   assert(_mesh);
 
   const int commRank = _mesh->commRank();
@@ -105,6 +111,8 @@
     _setMaterials(materialIds);
   }
   _distributeGroups();
+
+  PYLITH_METHOD_END;
 } // read
 
 // ----------------------------------------------------------------------
@@ -112,11 +120,15 @@
 void
 pylith::meshio::MeshIOCubit::_write(void) const
 { // write
+  PYLITH_METHOD_BEGIN;
+
   ExodusII exofile(_filename.c_str());
 
   _writeDimensions(exofile);
   _writeVariables(exofile);
   _writeAttributes(exofile);
+
+  PYLITH_METHOD_END;
 } // write
 
 // ----------------------------------------------------------------------
@@ -127,10 +139,12 @@
 					   int* numVertices, 
 					   int* numDims) const
 { // _readVertices
-  assert(0 != coordinates);
-  assert(0 != numVertices);
-  assert(0 != numDims);
+  PYLITH_METHOD_BEGIN;
 
+  assert(coordinates);
+  assert(numVertices);
+  assert(numDims);
+
   journal::info_t info("meshiocubit");
     
   // Space dimension
@@ -173,6 +187,8 @@
 	(*coordinates)[iVertex*(*numDims)+i] = buffer[iVertex];
     } // for
   } // else
+
+  PYLITH_METHOD_END;
 } // _readVertices
 
 // ----------------------------------------------------------------------
@@ -184,11 +200,13 @@
 					int* numCells, 
 					int* numCorners) const
 { // _readCells
-  assert(0 != cells);
-  assert(0 != materialIds);
-  assert(0 != numCells);
-  assert(0 != numCorners);
+  PYLITH_METHOD_BEGIN;
 
+  assert(cells);
+  assert(materialIds);
+  assert(numCells);
+  assert(numCorners);
+
   journal::info_t info("meshiocubit");
 
   *numCells = exofile.getDim("num_elem");
@@ -243,6 +261,8 @@
   } // for
 
   *cells -= 1; // use zero index
+
+  PYLITH_METHOD_END;
 } // _readCells
 
 // ----------------------------------------------------------------------
@@ -250,6 +270,8 @@
 void
 pylith::meshio::MeshIOCubit::_readGroups(ExodusII& exofile)
 { // _readGroups
+  PYLITH_METHOD_BEGIN;
+
   journal::info_t info("meshiocubit");
 
   const int numGroups = exofile.getDim("num_node_sets");
@@ -300,6 +322,8 @@
       _setGroup(name.str().c_str(), type, points);
     } // if/else
   } // for  
+
+  PYLITH_METHOD_END;
 } // _readGroups
 
 // ----------------------------------------------------------------------
@@ -334,7 +358,9 @@
 					  const int numCorners,
 					  const int meshDim)
 { // _orientCells
-  assert(0 != cells);
+  PYLITH_METHOD_BEGIN;
+
+  assert(cells);
   assert(cells->size() == numCells*numCorners);
 
   if (2 == meshDim && 4 == numCorners) { // QUAD4
@@ -429,6 +455,7 @@
     } // for
   } // if/else
 
+  PYLITH_METHOD_END;
 } // _orientCells
   
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOLagrit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOLagrit.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIOLagrit.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -55,7 +55,11 @@
 void
 pylith::meshio::MeshIOLagrit::deallocate(void)
 { // deallocate
+  PYLITH_METHOD_BEGIN;
+
   MeshIO::deallocate();
+
+  PYLITH_METHOD_END;
 } // deallocate
   
 // ----------------------------------------------------------------------
@@ -63,6 +67,8 @@
 void
 pylith::meshio::MeshIOLagrit::_read(void)
 { // _read
+  PYLITH_METHOD_BEGIN;
+
   const int commRank = _mesh->commRank();
   int meshDim = 0;
   int spaceDim = 0;
@@ -109,6 +115,8 @@
       _setGroup(groups[iGroup].name, type, groups[iGroup].points);
   }
   _distributeGroups();
+
+  PYLITH_METHOD_END;
 } // _read
 
 // ----------------------------------------------------------------------
@@ -128,7 +136,9 @@
 						const int numCorners,
 						const int meshDim)
 { // _orientCellsAscii
-  assert(0 != cells);
+  PYLITH_METHOD_BEGIN;
+
+  assert(cells);
   assert(cells->size() == numCells*numCorners);
 
   if (3 == meshDim && 4 == numCorners) // TET
@@ -139,6 +149,8 @@
       (*cells)[i1] = (*cells)[i2];
       (*cells)[i2] = tmp;
     } // for
+
+  PYLITH_METHOD_END;
 } // _orientCellsAscii
   
 // ----------------------------------------------------------------------
@@ -150,11 +162,15 @@
 						 const int numCorners,
 						 const int meshDim)
 { // _orientCellsBinary
-  assert(0 != cells);
+  PYLITH_METHOD_BEGIN;
+
+  assert(cells);
   assert(cells->size() == numCells*numCorners);
 
   if (3 == meshDim && 4 == numCorners)  // TET
     ; // do nothing
+
+  PYLITH_METHOD_END;
 } // _orientCellsBinary
   
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFile.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFile.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFile.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -22,6 +22,8 @@
 
 #include "PsetFileAscii.hh"
 
+#include "pylith/utils/petscerror.h" // USES PYLITH_METHOD_BEGIN/END
+
 #include <fstream> // uses std::fstream
 #include <sstream> // uses std::ostringstream
 #include <cstring> // uses strcmp()
@@ -42,6 +44,8 @@
 bool
 pylith::meshio::PsetFile::isAscii(const char* filename)
 { // isAscii
+  PYLITH_METHOD_BEGIN;
+
   std::ifstream fin(filename);
   if (!(fin.is_open() && fin.good())) {
     std::ostringstream msg;
@@ -52,7 +56,9 @@
   char buffer[headerLen];
   fin.get(buffer, headerLen, '\n');
   fin.close();
-  return (0 == strcmp(PsetFileAscii::header(), buffer)) ? true : false;
+
+  const bool result = (0 == strcmp(PsetFileAscii::header(), buffer)) ? true : false;
+  PYLITH_METHOD_RETURN(result);
 } // isAscii
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFileAscii.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFileAscii.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFileAscii.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -21,6 +21,7 @@
 #include "PsetFileAscii.hh" // implementation of class methods
 
 #include "pylith/utils/array.hh" // USES scalar_array, int_array
+#include "pylith/utils/petscerror.h" // USES PYLITH_METHOD_BEGIN/END
 
 #include "journal/info.h" // USES journal::info_t
 
@@ -52,8 +53,10 @@
 void
 pylith::meshio::PsetFileAscii::read(std::vector<Pset>* groups)
 { // read
-  assert(0 != groups);
+  PYLITH_METHOD_BEGIN;
 
+  assert(groups);
+
   journal::info_t info("psetfile");
 
   std::ifstream fin(_filename.c_str(), std::ios::in);
@@ -80,6 +83,8 @@
        << "Reading " << numGroups << " point sets from file." << journal::endl;
   for (int iGroup=0; iGroup < numGroups; ++iGroup)
     _readPset(fin, &(*groups)[iGroup]);
+
+  PYLITH_METHOD_END;
 } // read
 
 // ----------------------------------------------------------------------
@@ -87,6 +92,8 @@
 void
 pylith::meshio::PsetFileAscii::write(const std::vector<Pset>& groups)
 { // write
+  PYLITH_METHOD_BEGIN;
+
   journal::info_t info("psetfile");
 
   std::ofstream fout(_filename.c_str(), std::ios::out);
@@ -112,16 +119,20 @@
        << "Writing " << numGroups << " point sets to file." << journal::endl;
   for (int iGroup=0; iGroup < numGroups; ++iGroup)
     _writePset(fout, groups[iGroup]);
+
+  PYLITH_METHOD_END;
 } // write
 
 // ----------------------------------------------------------------------
 void
 pylith::meshio::PsetFileAscii::_readHeader(std::ifstream& fin)
 { // _readHeader
+  PYLITH_METHOD_BEGIN;
+
   const int headerLen = strlen(_HEADER)+1;
   char buffer[headerLen];
   fin.get(buffer, headerLen, '\n');
-  if (0 != strcmp(_HEADER, buffer)) {
+  if (strcmp(_HEADER, buffer)) {
     std::ostringstream msg;
     msg
       << "Header in ASCII Pset file '" << buffer
@@ -129,13 +140,19 @@
       << _HEADER << "'";
     throw std::runtime_error(msg.str());
   } // if
+
+  PYLITH_METHOD_END;
 } // _readHeader
 
 // ----------------------------------------------------------------------
 void
 pylith::meshio::PsetFileAscii::_writeHeader(std::ofstream& fout)
 { // _writeHeader
+  PYLITH_METHOD_BEGIN;
+
   fout << _HEADER << std::endl;
+
+  PYLITH_METHOD_END;
 } // _writeHeader
 
 // ----------------------------------------------------------------------
@@ -143,8 +160,10 @@
 pylith::meshio::PsetFileAscii::_readPset(std::ifstream& fin,
 					 Pset* group)
 { // _readPset
-  assert(0 != group);
+  PYLITH_METHOD_BEGIN;
 
+  assert(group);
+
   journal::info_t info("psetfile");
 
   int size = 0;
@@ -161,6 +180,8 @@
 
   info << journal::at(__HERE__)
        << "Done." << journal::endl;
+
+  PYLITH_METHOD_END;
 } // _readPset
 
 // ----------------------------------------------------------------------
@@ -168,6 +189,8 @@
 pylith::meshio::PsetFileAscii::_writePset(std::ofstream& fout,
 					  const Pset& group)
 { // _writePset
+  PYLITH_METHOD_BEGIN;
+
   journal::info_t info("psetfile");
   const int size = group.points.size();
 
@@ -188,6 +211,8 @@
 
   info << journal::at(__HERE__)
        << "Done." << journal::endl;
+
+  PYLITH_METHOD_END;
 } // _writePset
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFileBinary.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFileBinary.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/PsetFileBinary.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -23,6 +23,7 @@
 #include "BinaryIO.hh" // USES readString()
 
 #include "pylith/utils/array.hh" // USES scalar_array, int_array
+#include "pylith/utils/petscerror.h" // USES PYLITH_METHOD_BEGIN/END
 
 #include "journal/info.h" // USES journal::info_t
 
@@ -60,8 +61,10 @@
 void
 pylith::meshio::PsetFileBinary::read(std::vector<Pset>* groups)
 { // read
-  assert(0 != groups);
+  PYLITH_METHOD_BEGIN;
 
+  assert(groups);
+
   journal::info_t info("psetfile");
 
   std::ifstream fin(_filename.c_str(), std::ios::in);
@@ -109,6 +112,8 @@
     for (int iGroup=0; iGroup < numGroups; ++iGroup)
       _readPset64(fin, &(*groups)[iGroup]);
   } // if/else
+
+  PYLITH_METHOD_END;
 } // read
 
 // ----------------------------------------------------------------------
@@ -116,6 +121,8 @@
 void
 pylith::meshio::PsetFileBinary::write(const std::vector<Pset>& groups)
 { // write
+  PYLITH_METHOD_BEGIN;
+
   journal::info_t info("psetfile");
 
   std::ofstream fout(_filename.c_str(), std::ios::out);
@@ -157,12 +164,16 @@
     for (int iGroup=0; iGroup < numGroups; ++iGroup)
       _writePset64(fout, groups[iGroup]);
   } // if/else
+
+  PYLITH_METHOD_END;
 } // write
 
 // ----------------------------------------------------------------------
 void
 pylith::meshio::PsetFileBinary::_readHeader(std::ifstream& fin)
 { // _readHeader
+  PYLITH_METHOD_BEGIN;
+
   std::string extra = BinaryIO::readString(fin, _recordHeaderSize);
 
   std::string header = BinaryIO::readString(fin, strlen(_HEADER));
@@ -175,13 +186,19 @@
       << "' does not match anticipated header '" << headerE << "'.";
     throw std::runtime_error(msg.str());
   } // if
+
+  PYLITH_METHOD_END;
 } // _readHeader
 
 // ----------------------------------------------------------------------
 void
 pylith::meshio::PsetFileBinary::_writeHeader(std::ofstream& fout)
 { // _writeHeader
+  PYLITH_METHOD_BEGIN;
+
   fout.write((char*) _HEADER, strlen(_HEADER));
+
+  PYLITH_METHOD_END;
 } // _writeHeader
 
 // ----------------------------------------------------------------------
@@ -189,8 +206,10 @@
 pylith::meshio::PsetFileBinary::_readPset32(std::ifstream& fin,
 					    Pset* group)
 { // _readPset32
-  assert(0 != group);
+  PYLITH_METHOD_BEGIN;
 
+  assert(group);
+
   journal::info_t info("psetfile");
 
   group->name = BinaryIO::readString(fin, 32);
@@ -220,6 +239,8 @@
 
   info << journal::at(__HERE__)
        << "Done." << journal::endl;
+
+  PYLITH_METHOD_END;
 } // _readPset32
 
 // ----------------------------------------------------------------------
@@ -227,6 +248,8 @@
 pylith::meshio::PsetFileBinary::_writePset32(std::ofstream& fout,
 					   const Pset& group)
 { // _writePset32
+  PYLITH_METHOD_BEGIN;
+
   journal::info_t info("psetfile");
   const int size = group.points.size();
   info << journal::at(__HERE__)
@@ -255,6 +278,8 @@
 
   info << journal::at(__HERE__)
        << "Done." << journal::endl;
+
+  PYLITH_METHOD_END;
 } // _writePset32
 
 // ----------------------------------------------------------------------
@@ -262,8 +287,10 @@
 pylith::meshio::PsetFileBinary::_readPset64(std::ifstream& fin,
 					    Pset* group)
 { // _readPset64
-  assert(0 != group);
+  PYLITH_METHOD_BEGIN;
 
+  assert(group);
+
   journal::info_t info("psetfile");
 
   group->name = BinaryIO::readString(fin, 32);
@@ -297,6 +324,8 @@
 
   info << journal::at(__HERE__)
        << "Done." << journal::endl;
+
+  PYLITH_METHOD_END;
 } // _readPset64
 
 // ----------------------------------------------------------------------
@@ -304,6 +333,8 @@
 pylith::meshio::PsetFileBinary::_writePset64(std::ofstream& fout,
 					     const Pset& group)
 { // _writePset64
+  PYLITH_METHOD_BEGIN;
+
   journal::info_t info("psetfile");
   const int size = group.points.size();
   info << journal::at(__HERE__)
@@ -330,6 +361,8 @@
 
   info << journal::at(__HERE__)
        << "Done." << journal::endl;
+
+  PYLITH_METHOD_END;
 } // _writePset64
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/Xdmf.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/Xdmf.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/Xdmf.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -23,7 +23,8 @@
 #include "HDF5.hh" // USES HDF5
 
 #include "pylith/topology/FieldBase.hh" // USES FieldBase enums
-#include "pylith/utils/array.hh" // USES PylithScalar_arra
+#include "pylith/utils/array.hh" // USES PylithScalar_array
+#include "pylith/utils/petscerror.h" // USES PYLITH_METHOD_BEGIN/END
 
 #include <string> // USES std::string
 #include <stdexcept> // USES std::runtime_error
@@ -49,6 +50,8 @@
 pylith::meshio::Xdmf::write(const char* filenameXdmf,
 			    const char* filenameHDF5)
 { // write
+  PYLITH_METHOD_BEGIN;
+
   assert(filenameXdmf);
   assert(filenameHDF5);
 
@@ -212,6 +215,8 @@
     << "  </Domain>\n"
     << "</Xdmf>\n";
   _file.close();
+
+  PYLITH_METHOD_END;
 } // write
 
 // ----------------------------------------------------------------------
@@ -220,6 +225,8 @@
 pylith::meshio::Xdmf::_getTimeStamps(scalar_array* timeStamps,
 				     HDF5& h5)
 { // _getTimeStamps
+  PYLITH_METHOD_BEGIN;
+
   assert(timeStamps);
 
   hsize_t* dims = 0;
@@ -250,6 +257,8 @@
 
   delete[] dims; dims = 0;
   delete[] t; t = 0;
+
+  PYLITH_METHOD_END;
 } // _getTimeStamps
 
 // ----------------------------------------------------------------------
@@ -258,6 +267,8 @@
 pylith::meshio::Xdmf::_getFieldMetadata(std::vector<FieldMetadata>* metadata,
 					HDF5& h5)
 { // _getFieldMetadata
+  PYLITH_METHOD_BEGIN;
+
   assert(metadata);
 
   string_vector fieldNames;
@@ -364,6 +375,8 @@
 	      << std::endl;
   } // for
 #endif
+
+  PYLITH_METHOD_END;
 } // _getFieldMetadata
 
 // ----------------------------------------------------------------------
@@ -372,6 +385,8 @@
 pylith::meshio::Xdmf::_writeDomainCells(const int numCells,
 					const int numCorners)
 { // _writeDomainCells
+  PYLITH_METHOD_BEGIN;
+
   assert(_file.is_open() && _file.good());
   _file
     << "    <DataItem Name=\"cells\"\n"
@@ -382,6 +397,7 @@
     << "      &HeavyData;:/topology/cells\n"
     << "    </DataItem>\n";
 
+  PYLITH_METHOD_END;
 } // _writeDomainCells
 
 // ----------------------------------------------------------------------
@@ -390,6 +406,8 @@
 pylith::meshio::Xdmf::_writeDomainVertices(const int numVertices,
 					   const int spaceDim)
 { // _writeDomainVertices
+  PYLITH_METHOD_BEGIN;
+
   assert(_file.is_open() && _file.good());
 
   _file
@@ -398,6 +416,8 @@
     << "	      Dimensions=\"" << numVertices << " " << spaceDim << "\">\n"
     << "      &HeavyData;:/geometry/vertices\n"
     << "    </DataItem>\n";
+
+  PYLITH_METHOD_END;
 } // _writeDomainVertices
 
 // ----------------------------------------------------------------------
@@ -405,6 +425,8 @@
 void
 pylith::meshio::Xdmf::_writeTimeStamps(const scalar_array& timeStamps)
 { // _writeTimeStamps
+  PYLITH_METHOD_BEGIN;
+
   assert(_file.is_open() && _file.good());
 
   const int numTimeStamps = timeStamps.size();
@@ -419,6 +441,8 @@
     << "\n"
     << "        </DataItem>\n"
     << "      </Time>\n";
+
+  PYLITH_METHOD_END;
 } // _writeTimeStamps
 
 // ----------------------------------------------------------------------
@@ -427,6 +451,8 @@
 pylith::meshio::Xdmf::_writeGridTopology(const char* cellType,
 					 const int numCells)
 { // _writeGridTopology
+  PYLITH_METHOD_BEGIN;
+
   assert(_file.is_open() && _file.good());
 
   _file
@@ -437,6 +463,8 @@
     << "	    /Xdmf/Domain/DataItem[@Name=\"cells\"]\n"
     << "	  </DataItem>\n"
     << "	</Topology>\n";
+
+  PYLITH_METHOD_END;
 } // _writeGridTopology
 
 // ----------------------------------------------------------------------
@@ -444,6 +472,8 @@
 void
 pylith::meshio::Xdmf::_writeGridGeometry(const int spaceDim)
 { // _writeGridGeometry
+  PYLITH_METHOD_BEGIN;
+
   assert(_file.is_open() && _file.good());
   assert(2 == spaceDim || 3 == spaceDim);
 
@@ -455,6 +485,8 @@
     << "	    /Xdmf/Domain/DataItem[@Name=\"vertices\"]\n"
     << "	  </DataItem>\n"
     << "	</Geometry>\n";
+
+  PYLITH_METHOD_END;
 } // _writeGridGeometry
 
 // ----------------------------------------------------------------------
@@ -463,6 +495,8 @@
 pylith::meshio::Xdmf::_writeGridAttribute(const FieldMetadata& metadata,
 					  const int iTime)
 { // _writeGridAttribute
+  PYLITH_METHOD_BEGIN;
+
   assert(_file.is_open() && _file.good());
 
   std::string h5FullName = "";
@@ -507,6 +541,8 @@
     << "	    </DataItem>\n"
     << "	  </DataItem>\n"
     << "	</Attribute>\n";
+
+  PYLITH_METHOD_END;
 } // _writeGridAttribute
 
 // ----------------------------------------------------------------------
@@ -517,6 +553,8 @@
 						   const int component,
 						   const int spaceDim)
 { // _writeGridAttribute
+  PYLITH_METHOD_BEGIN;
+
   assert(_file.is_open() && _file.good());
 
   std::string h5FullName = "";
@@ -659,6 +697,8 @@
     << "	    </DataItem>\n"
     << "	  </DataItem>\n"
     << "	</Attribute>\n";
+
+  PYLITH_METHOD_END;
 } // _writeGridAttribute
 
 

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -21,6 +21,8 @@
 #include "TestMeshIO.hh" // Implementation of class methods
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
 #include "pylith/meshio/MeshIO.hh" // USES MeshIO
 
 #include "pylith/utils/array.hh" // USES int_array
@@ -34,18 +36,30 @@
 #include <stdexcept> // USES std::logic_error
 
 // ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
+// Setup testing data.
+void
+pylith::meshio::TestMeshIO::setUp(void)
+{ // setUp
+  _mesh = 0;
+} // setUp
 
 // ----------------------------------------------------------------------
+// Tear down testing data.
+void
+pylith::meshio::TestMeshIO::tearDown(void)
+{ // tearDown
+  delete _mesh; _mesh = 0;
+} // tearDown
+
+// ----------------------------------------------------------------------
 // Get simple mesh for testing I/O.
-pylith::topology::Mesh*
+void
 pylith::meshio::TestMeshIO::_createMesh(const MeshData& data)
 { // _createMesh
   PYLITH_METHOD_BEGIN;
 
   // buildTopology() requires zero based index
-  CPPUNIT_ASSERT(true == data.useIndexZero);
+  CPPUNIT_ASSERT_EQUAL(true, data.useIndexZero);
 
   CPPUNIT_ASSERT(data.vertices);
   CPPUNIT_ASSERT(data.cells);
@@ -57,10 +71,11 @@
     CPPUNIT_ASSERT(data.groupTypes);
   } // if
 
-  topology::Mesh* mesh = new topology::Mesh(data.cellDim);CPPUNIT_ASSERT(mesh);
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();CPPUNIT_ASSERT(!sieveMesh.isNull());
-  ALE::Obj<SieveMesh::sieve_type> sieve = new SieveMesh::sieve_type(mesh->comm());CPPUNIT_ASSERT(!sieve.isNull());
+  delete _mesh; _mesh = new topology::Mesh(data.cellDim);CPPUNIT_ASSERT(_mesh);
 
+  const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();CPPUNIT_ASSERT(!sieveMesh.isNull());
+  ALE::Obj<SieveMesh::sieve_type> sieve = new SieveMesh::sieve_type(_mesh->comm());CPPUNIT_ASSERT(!sieve.isNull());
+
   // Cells and vertices
   const bool interpolate = false;
   ALE::Obj<SieveFlexMesh::sieve_type> s = new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());CPPUNIT_ASSERT(s);
@@ -75,8 +90,8 @@
   PetscBool interpolateMesh = interpolate ? PETSC_TRUE : PETSC_FALSE;
   PetscErrorCode err;
 
-  err = DMPlexCreateFromCellList(mesh->comm(), data.cellDim, data.numCells, data.numVertices, data.numCorners, interpolateMesh, data.cells, data.spaceDim, data.vertices, &dmMesh);CHECK_PETSC_ERROR(err);
-  mesh->setDMMesh(dmMesh);
+  err = DMPlexCreateFromCellList(_mesh->comm(), data.cellDim, data.numCells, data.numVertices, data.numCorners, interpolateMesh, data.cells, data.spaceDim, data.vertices, &dmMesh);CHECK_PETSC_ERROR(err);
+  _mesh->setDMMesh(dmMesh);
 
   // Material ids
   const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
@@ -126,199 +141,89 @@
   spatialdata::units::Nondimensional normalizer;
   cs.setSpaceDim(data.spaceDim);
   cs.initialize();
-  mesh->coordsys(&cs);
-  mesh->nondimensionalize(normalizer);
+  _mesh->coordsys(&cs);
+  _mesh->nondimensionalize(normalizer);
 
-  PYLITH_METHOD_RETURN(mesh);
+  PYLITH_METHOD_END;
 } // _createMesh
 
 // ----------------------------------------------------------------------
 // Check values in mesh against data.
 void
-pylith::meshio::TestMeshIO::_checkVals(const topology::Mesh& mesh,
-				       const MeshData& data)
+pylith::meshio::TestMeshIO::_checkVals(const MeshData& data)
 { // _checkVals
   PYLITH_METHOD_BEGIN;
 
-  typedef SieveMesh::int_section_type::chart_type chart_type;
+  CPPUNIT_ASSERT(_mesh);
 
   // Check mesh dimension
-  CPPUNIT_ASSERT_EQUAL(data.cellDim, mesh.dimension());
+  CPPUNIT_ASSERT_EQUAL(data.cellDim, _mesh->dimension());
+  const int spaceDim = data.spaceDim;
 
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();CPPUNIT_ASSERT(!sieveMesh.isNull());
+  PetscDM dmMesh = _mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
 
-  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
-
   // Check vertices
-  const ALE::Obj<SieveMesh::label_sequence>& vertices = 
-    sieveMesh->depthStratum(0);
-  const ALE::Obj<RealSection>& coordsField =
-    sieveMesh->getRealSection("coordinates");
-  const int numVertices = vertices->size();
-  CPPUNIT_ASSERT(!vertices.isNull());
-  CPPUNIT_ASSERT(!coordsField.isNull());
-  CPPUNIT_ASSERT_EQUAL(data.numVertices, numVertices);
-  CPPUNIT_ASSERT_EQUAL(data.spaceDim, coordsField->getFiberDimension(*vertices->begin()));
-  int i = 0;
-  const int spaceDim = data.spaceDim;
-  for(SieveMesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
-    const PylithScalar* vertexCoords = coordsField->restrictPoint(*v_iter);CPPUNIT_ASSERT(vertexCoords);
-    const PylithScalar tolerance = 1.0e-06;
-    for (int iDim=0; iDim < spaceDim; ++iDim)
-      if (data.vertices[i] < 1.0) {
-        CPPUNIT_ASSERT_DOUBLES_EQUAL(data.vertices[i++], vertexCoords[iDim],
-				     tolerance);
-      } else {
-        CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vertexCoords[iDim]/data.vertices[i++],
-				     tolerance);
-      } // if/else
-  } // for
-  PetscSection coordSection = NULL;
-  PetscVec coordVec = NULL;
-  PetscScalar *coords = NULL;
-  PetscInt vStart, vEnd;
-  PetscErrorCode err;
+  topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
 
-  err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);CPPUNIT_ASSERT(coordSection);
-  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);CPPUNIT_ASSERT(coordVec);
-  err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    PetscInt dof, off;
+  CPPUNIT_ASSERT_EQUAL(data.numVertices, verticesStratum.size());
 
-    err = PetscSectionGetDof(coordSection, v, &dof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(coordSection, v, &off);CHECK_PETSC_ERROR(err);
-    assert(spaceDim == dof);
-    const PylithScalar tolerance = 1.0e-06;
-    for(PetscInt d = 0; d < spaceDim; ++d) {
-      const PetscInt eoff = (v-vStart)*spaceDim;
-      if (data.vertices[eoff+d] < 1.0) {
-        CPPUNIT_ASSERT_DOUBLES_EQUAL(data.vertices[eoff+d], coords[off+d], tolerance);
+  topology::CoordsVisitor coordsVisitor(dmMesh);
+  const PetscScalar* coordsArray = coordsVisitor.localArray();
+  const PylithScalar tolerance = 1.0e-06;
+  for (PetscInt v = vStart, index = 0; v < vEnd; ++v) {
+    const PetscInt off = coordsVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, coordsVisitor.sectionDof(v));
+
+    for (int iDim=0; iDim < spaceDim; ++iDim, ++index) {
+      if (data.vertices[index] < 1.0) {
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(data.vertices[index], coordsArray[off+iDim], tolerance);
       } else {
-        CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, coords[off+d]/data.vertices[eoff+d], tolerance);
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, coordsArray[off+iDim]/data.vertices[index], tolerance);
       } // if/else
     } // for
   } // for
-  err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
 
-  // check cells
-  const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
-  const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
+  // Check cells
+  topology::Stratum cellsStratum(dmMesh, topology::Stratum::HEIGHT, 0);
+  const PetscInt cStart = cellsStratum.begin();
+  const PetscInt cEnd = cellsStratum.end();
+  const PetscInt numCells = cellsStratum.size();
 
-  const int numCells = cells->size();
   CPPUNIT_ASSERT_EQUAL(data.numCells, numCells);
-  const int numCorners = sieveMesh->getNumCellCorners();
-  CPPUNIT_ASSERT_EQUAL(data.numCorners, numCorners);
-
-  ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> pV(sieve->getMaxConeSize());
   const int offset = numCells;
-  i = 0;
-  for(SieveMesh::label_sequence::iterator e_iter = cells->begin();
-      e_iter != cells->end();
-      ++e_iter) {
-    sieve->cone(*e_iter, pV);
-    const SieveMesh::point_type *cone = pV.getPoints();
-    for(int p = 0; p < pV.getSize(); ++p, ++i) {
-      CPPUNIT_ASSERT_EQUAL(data.cells[i], cone[p]-offset);
-    } // for
-    pV.clear();
-  } // for
-  PetscInt cStart, cEnd;
-
-  err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT_EQUAL(data.numCells, cEnd-cStart);
-  for(PetscInt c = cStart; c < cEnd; ++c) {
-    const PetscInt *cone;
-    PetscInt        coneSize;
-
+  const PetscInt* cone = NULL;
+  PetscInt coneSize = 0;
+  PetscErrorCode err = 0;
+  for(PetscInt c = cStart, index = 0; c < cEnd; ++c) {
     err = DMPlexGetConeSize(dmMesh, c, &coneSize);CHECK_PETSC_ERROR(err);
     err = DMPlexGetCone(dmMesh, c, &cone);CHECK_PETSC_ERROR(err);
     CPPUNIT_ASSERT_EQUAL(data.numCorners, coneSize);
-    for(PetscInt p = 0; p < coneSize; ++p) {
-      CPPUNIT_ASSERT_EQUAL(data.cells[(c-cStart)*coneSize+p], cone[p]-offset);
+    for(PetscInt p = 0; p < coneSize; ++p, ++index) {
+      CPPUNIT_ASSERT_EQUAL(data.cells[index], cone[p]-offset);
     } // for
   } // for
 
   // check materials
-  const ALE::Obj<SieveMesh::label_type>& labelMaterials = sieveMesh->getLabel("material-id");
-  const int idDefault = -999;
-  const int size = numCells;
-  int_array materialIds(size);
-  i = 0;
-  for(SieveMesh::label_sequence::iterator e_iter = cells->begin();
-      e_iter != cells->end();
-      ++e_iter)
-    materialIds[i++] = sieveMesh->getValue(labelMaterials, *e_iter, idDefault);
-  
-  for (int iCell=0; iCell < numCells; ++iCell)
-    CPPUNIT_ASSERT_EQUAL(data.materialIds[iCell], materialIds[iCell]);
-
+  PetscInt matId = 0;
   for(PetscInt c = cStart; c < cEnd; ++c) {
-    PetscInt matId;
-
     err = DMPlexGetLabelValue(dmMesh, "material-id", c, &matId);CHECK_PETSC_ERROR(err);
     CPPUNIT_ASSERT_EQUAL(data.materialIds[c-cStart], matId);
   } // for
 
   // Check groups
-#if 0 // SIEVE STUFF
-  const ALE::Obj<std::set<std::string> >& groupNames = 
-    sieveMesh->getIntSections();
-  if (data.numGroups > 0) {
-    CPPUNIT_ASSERT(!groupNames.isNull());
-    CPPUNIT_ASSERT_EQUAL(data.numGroups, int(groupNames->size()));
-  } // if
-  int iGroup = 0;
-  int index = 0;
-  for (std::set<std::string>::const_iterator name=groupNames->begin();
-       name != groupNames->end();
-       ++name, ++iGroup) {
-    const ALE::Obj<SieveMesh::int_section_type>& groupField = 
-      sieveMesh->getIntSection(*name);
-    CPPUNIT_ASSERT(!groupField.isNull());
-    const chart_type& chart = groupField->getChart();
-    SieveMesh::point_type firstPoint;
-    for(chart_type::const_iterator c_iter = chart.begin();
-	c_iter != chart.end();
-	++c_iter) {
-      if (groupField->getFiberDimension(*c_iter)) {
-	firstPoint = *c_iter;
-	break;
-      } // if
-    } // for
-    std::string groupType = 
-      (sieveMesh->height(firstPoint) == 0) ? "cell" : "vertex";
-    const int numPoints = groupField->size();
-    int_array points(numPoints);
-    int i = 0;
-    const int offset = ("vertex" == groupType) ? numCells : 0;
-    for(chart_type::const_iterator c_iter = chart.begin();
-	c_iter != chart.end();
-	++c_iter)
-      if (groupField->getFiberDimension(*c_iter))
-	points[i++] = *c_iter - offset;
-    
-    CPPUNIT_ASSERT_EQUAL(std::string(data.groupNames[iGroup]), *name);
-    CPPUNIT_ASSERT_EQUAL(std::string(data.groupTypes[iGroup]), groupType);
-    CPPUNIT_ASSERT_EQUAL(data.groupSizes[iGroup], numPoints);
-    for (int i=0; i < numPoints; ++i)
-      CPPUNIT_ASSERT_EQUAL(data.groups[index++], points[i]);
-  } // for
-#endif
-
-  PetscInt numLabels, pStart, pEnd;
-
+  PetscInt numGroups, pStart, pEnd;
   err = DMPlexGetChart(dmMesh, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetNumLabels(dmMesh, &numLabels);CHECK_PETSC_ERROR(err);
-  numLabels -= 2; /* Remove depth and material labels */
-  CPPUNIT_ASSERT_EQUAL(data.numGroups, numLabels);
-  PetscInt iGroup = 0;
+  err = DMPlexGetNumLabels(dmMesh, &numGroups);CHECK_PETSC_ERROR(err);
+  numGroups -= 2; // Remove depth and material labels.
+  CPPUNIT_ASSERT_EQUAL(data.numGroups, numGroups);
   PetscInt index  = 0;
-  for(PetscInt l = numLabels-1; l >= 0; --l, ++iGroup) {
-    const char *name;
-    PetscInt    firstPoint;
+  for(PetscInt iGroup = 0, iLabel = numGroups-1; iGroup < numGroups; ++iGroup, --iLabel) {
+    const char *name = NULL;
+    PetscInt firstPoint = 0;
 
-    err = DMPlexGetLabelName(dmMesh, l, &name);CHECK_PETSC_ERROR(err);
+    err = DMPlexGetLabelName(dmMesh, iLabel, &name);CHECK_PETSC_ERROR(err);
     CPPUNIT_ASSERT_EQUAL(std::string(data.groupNames[iGroup]), std::string(name));
     for(PetscInt p = pStart; p < pEnd; ++p) {
       PetscInt val;

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.hh	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.hh	2013-04-13 16:36:29 UTC (rev 21859)
@@ -47,6 +47,15 @@
 class pylith::meshio::TestMeshIO : public CppUnit::TestFixture
 { // class TestMeshIO
 
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public :
+  
+  /// Setup testing data.
+  void setUp(void);
+
+  /// Tear down testing data.
+  void tearDown(void);
+
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :
 
@@ -56,15 +65,13 @@
    *
    * @returns PyLith mesh
    */
-  topology::Mesh* _createMesh(const MeshData& data);
+  void _createMesh(const MeshData& data);
 
   /** Check values in mesh against data.
    *
-   * @param mesh PyLith mesh
    * @param data Mesh data
    */
-  void _checkVals(const topology::Mesh& mesh,
-		  const MeshData& data);
+  void _checkVals(const MeshData& data);
 
   /** Test debug().
    *
@@ -78,6 +85,11 @@
    */
   void _testInterpolate(MeshIO& iohandler);
 
+// PROTECTED MEMEBERS ////////////////////////////////////////////////////
+protected :
+
+  topology::Mesh* _mesh; ///< Finite-element mesh.
+
 }; // class TestMeshIO
 
 #endif // pylith_meshio_testmeshio_hh

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOAscii.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOAscii.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOAscii.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -211,20 +211,19 @@
 { // _testWriteRead
   PYLITH_METHOD_BEGIN;
 
-  topology::Mesh* meshOut = _createMesh(data);
+  _createMesh(data);
 
   // Write mesh
   MeshIOAscii iohandler;
   iohandler.filename(filename);
-  iohandler.write(meshOut);
-  delete meshOut; meshOut = 0;
+  iohandler.write(_mesh);
 
   // Read mesh
-  topology::Mesh meshIn;
-  iohandler.read(&meshIn);
+  delete _mesh; _mesh = new topology::Mesh;
+  iohandler.read(_mesh);
 
   // Make sure meshIn matches data
-  _checkVals(meshIn, data);
+  _checkVals(data);
 
   PYLITH_METHOD_END;
 } // _testWriteRead
@@ -238,13 +237,13 @@
   PYLITH_METHOD_BEGIN;
 
   // Read mesh
-  topology::Mesh mesh;
   MeshIOAscii iohandler;
   iohandler.filename(filename);
-  iohandler.read(&mesh);
+  delete _mesh; _mesh = new topology::Mesh;
+  iohandler.read(_mesh);
 
   // Make sure mesh matches data
-  _checkVals(mesh, data);
+  _checkVals(data);
 
   PYLITH_METHOD_END;
 } // _testRead

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOCubit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOCubit.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOCubit.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -127,11 +127,11 @@
   iohandler.useNodesetNames(true);
 
   // Read mesh
-  topology::Mesh mesh;
-  iohandler.read(&mesh);
+  delete _mesh; _mesh = new topology::Mesh;
+  iohandler.read(_mesh);
 
-  // Make sure meshIn matches data
-  _checkVals(mesh, data);
+  // Make sure mesh matches data
+  _checkVals(data);
 } // _testRead
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOLagrit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOLagrit.cc	2013-04-13 15:47:19 UTC (rev 21858)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIOLagrit.cc	2013-04-13 16:36:29 UTC (rev 21859)
@@ -155,11 +155,11 @@
 #endif
 
   // Read mesh
-  topology::Mesh mesh;
-  iohandler.read(&mesh);
+  delete _mesh; _mesh = new topology::Mesh;
+  iohandler.read(_mesh);
 
-  // Make sure meshIn matches data
-  _checkVals(mesh, data);
+  // Make sure mesh matches data
+  _checkVals(data);
 } // _testRead
 
 // ----------------------------------------------------------------------



More information about the CIG-COMMITS mailing list