[cig-commits] r13936 - in short/3D/PyLith/branches/pylith-swig: libsrc/faults libsrc/meshio unittests/libtests/meshio

brad at geodynamics.org brad at geodynamics.org
Fri Jan 23 16:22:09 PST 2009


Author: brad
Date: 2009-01-23 16:22:09 -0800 (Fri, 23 Jan 2009)
New Revision: 13936

Modified:
   short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/meshio/Makefile.am
   short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIO.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIO.hh
   short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOAscii.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOAscii.hh
   short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOCubit.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOCubit.hh
   short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOLagrit.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOLagrit.hh
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/Makefile.am
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestCellFilterAvg.cc
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestCellFilterAvg.hh
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIO.cc
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIO.hh
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOAscii.cc
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOAscii.hh
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOCubit.cc
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOCubit.hh
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOLagrit.cc
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOLagrit.hh
Log:
Updated MeshIO to use PyLith Mesh object. Temporarily disabled tests of other meshio objects.

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -16,7 +16,9 @@
 
 #include "CohesiveTopology.hh" // USES CohesiveTopology::create()
 
-#include "pylith/meshio/MeshIOLagrit.hh" // USES MeshIOLagrit::readFault()
+#if defined(NEWPYLITHMEHS)
+#include "pylith/meshio/UCDFaultFile.hh" // USES UCDFaultFile::read()
+#endif
 
 #include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
 #include "pylith/utils/array.hh" // USES double_array
@@ -71,8 +73,10 @@
 
     //MPI_Bcast(&faultDim, 1, MPI_INT, 0, comm);
     faultMesh = new SubMesh(mesh->comm(), faultDim, mesh->debug());
-    pylith::meshio::MeshIOLagrit::readFault(_faultMeshFilename, mesh, 
-					    faultMesh, faultBd);
+#if defined(NEWPYLITHMESH)
+    meshio::UCDFaultFile::readFault(_faultMeshFilename, mesh, 
+				    faultMesh, faultBd);
+#endif
 
     // Get group of vertices associated with fault
     const ALE::Obj<int_section_type>& groupField = 

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/Makefile.am	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/Makefile.am	2009-01-24 00:22:09 UTC (rev 13936)
@@ -19,6 +19,7 @@
 	DataWriter.hh \
 	DataWriterVTK.hh \
 	DataWriterVTK.icc \
+	MeshBuilder.hh \
 	MeshIO.hh \
 	MeshIO.icc \
 	MeshIOAscii.hh \
@@ -28,7 +29,8 @@
 	OutputManager.hh \
 	OutputSolnSubset.hh \
 	VertexFilter.hh \
-	VertexFilterVecNorm.hh
+	VertexFilterVecNorm.hh \
+	UCDFaultFile.cc
 
 if ENABLE_CUBIT
   subpkginclude_HEADERS += \

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIO.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIO.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIO.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -14,6 +14,7 @@
 
 #include "MeshIO.hh" // implementation of class methods
 
+#include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/utils/array.hh" // USES double_array, int_array
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -46,39 +47,38 @@
 int
 pylith::meshio::MeshIO::getMeshDim(void) const
 { // getMeshDim
-  assert(0 != _mesh);
-  assert(!_mesh->isNull());
-  return (*_mesh)->getDimension();
+  return (0 != _mesh) ? _mesh->dimension() : 0;
 } // getMeshDim
 
 // ----------------------------------------------------------------------
+// Set scales used to nondimensionalize mesh.
+void
+pylith::meshio::MeshIO::normalizer(
+			       const spatialdata::units::Nondimensional& dim)
+{ // normalizer
+  if (0 == _normalizer)
+    _normalizer = new spatialdata::units::Nondimensional(dim);
+  else
+    *_normalizer = dim;
+} // normalizer
+
+// ----------------------------------------------------------------------
 // Read mesh from file.
 void 
-pylith::meshio::MeshIO::read(ALE::Obj<Mesh>* mesh)
+pylith::meshio::MeshIO::read(topology::Mesh* mesh)
 { // read
   assert(0 == _mesh);
 
   _mesh = mesh;
+  _mesh->debug(_debug);
   _read();
   _mesh = 0;
-  (*mesh)->getFactory()->clear();
 } // read
 
 // ----------------------------------------------------------------------
-// Set scales used to nondimensionalize mesh.
-void
-pylith::meshio::MeshIO::normalizer(const spatialdata::units::Nondimensional& dim)
-{ // normalizer
-  if (0 == _normalizer)
-    _normalizer = new spatialdata::units::Nondimensional(dim);
-  else
-    *_normalizer = dim;
-} // normalizer
-
-// ----------------------------------------------------------------------
 // Write mesh to file.
 void 
-pylith::meshio::MeshIO::write(ALE::Obj<Mesh>* mesh)
+pylith::meshio::MeshIO::write(topology::Mesh* const mesh)
 { // write
   assert(0 == _mesh);
 
@@ -88,261 +88,6 @@
 } // write
 
 // ----------------------------------------------------------------------
-// Set vertices and cells in mesh.
-void
-pylith::meshio::MeshIO::_buildMesh(double_array* coordinates,
-				   const int numVertices,
-				   const int spaceDim,
-				   const int_array& cells,
-				   const int numCells,
-				   const int numCorners,
-				   const int meshDim)
-{ // _buildMesh
-  assert(0 != _mesh);
-  assert(0 != coordinates);
-  MPI_Comm comm = PETSC_COMM_WORLD;
-  int      dim  = meshDim;
-  int      rank;
-
-  { // Check to make sure every vertex is in at least one cell.
-    // This is required by Sieve
-    std::vector<bool> vertexInCell(numVertices, false);
-    const int size = cells.size();
-    for (int i=0; i < size; ++i)
-      vertexInCell[cells[i]] = true;
-    int count = 0;
-    for (int i=0; i < numVertices; ++i)
-      if (!vertexInCell[i])
-	++count;
-    if (count > 0) {
-      std::ostringstream msg;
-      msg << "Mesh contains " << count
-	  << " vertices that are not in any cells.";
-      throw std::runtime_error(msg.str());
-    } // if
-  } // check
-
-  MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
-  // :BUG: This causes a memory leak.
-  *_mesh = new Mesh(PETSC_COMM_WORLD, dim);
-
-  assert(!_mesh->isNull());
-  (*_mesh)->setDebug(_debug);
-
-  ALE::Obj<sieve_type> sieve = new sieve_type((*_mesh)->comm());
-  (*_mesh)->setSieve(sieve);
-
-  MPI_Comm_rank(comm, &rank);
-  // Memory debugging
-  ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
-  logger.setDebug(_debug/2);
-
-  logger.stagePush("MeshCreation");
-  if (!rank) {
-    assert(coordinates->size() == numVertices*spaceDim);
-    assert(cells.size() == numCells*numCorners);
-    if (!_interpolate) {
-      // Create the ISieve
-      sieve->setChart(Mesh::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);
-      }
-      delete [] cone;
-      // Symmetrize to fill up supports
-      sieve->symmetrize();
-    } else {
-      // Same old thing
-      ALE::Obj<ALE::Mesh::sieve_type> s = new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
-
-      ALE::SieveBuilder<ALE::Mesh>::buildTopology(s, meshDim, 
-                                                  numCells, 
-                                                  const_cast<int*>(&cells[0]), 
-                                                  numVertices, 
-                                                  _interpolate,
-                                                  numCorners);
-      std::map<Mesh::point_type,Mesh::point_type> renumbering;
-      ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
-    }
-    logger.stagePop();
-    logger.stagePush("MeshStratification");
-    if (!_interpolate) {
-      // Optimized stratification
-      const ALE::Obj<Mesh::label_type>& height = (*_mesh)->createLabel("height");
-      const ALE::Obj<Mesh::label_type>& depth  = (*_mesh)->createLabel("depth");
-
-#ifdef IMESH_NEW_LABELS
-      height->setChart((*_mesh)->getSieve()->getChart());
-      depth->setChart((*_mesh)->getSieve()->getChart());
-      for(int c = 0; c < numCells+numVertices; ++c) {
-        height->setConeSize(c, 1);
-        depth->setConeSize(c, 1);
-      }
-      if (numCells+numVertices) {height->setSupportSize(0, numCells+numVertices);}
-      if (numCells+numVertices) {depth->setSupportSize(0, numCells+numVertices);}
-      height->allocate();
-      depth->allocate();
-#endif
-      for(int c = 0; c < numCells; ++c) {
-        height->setCone(0, c);
-        depth->setCone(1, c);
-      }
-      for(int v = numCells; v < numCells+numVertices; ++v) {
-        height->setCone(1, v);
-        depth->setCone(0, v);
-      }
-#ifdef IMESH_NEW_LABELS
-      height->recalculateLabel();
-      depth->recalculateLabel();
-#endif
-      (*_mesh)->setHeight(1);
-      (*_mesh)->setDepth(1);
-    } else {
-      (*_mesh)->stratify();
-    }
-    logger.stagePop();
-  } else {
-    logger.stagePush("MeshStratification");
-    (*_mesh)->getSieve()->setChart(sieve_type::chart_type());
-    (*_mesh)->getSieve()->allocate();
-    (*_mesh)->stratify();
-    logger.stagePop();
-  }
-
-#if defined(ALE_MEM_LOGGING)
-  std::cout
-    << std::endl
-    << "MeshCreation " << logger.getNumAllocations("MeshCreation")
-    << " allocations " << logger.getAllocationTotal("MeshCreation")
-    << " bytes" << std::endl
-    
-    << "MeshCreation " << logger.getNumDeallocations("MeshCreation")
-    << " deallocations " << logger.getDeallocationTotal("MeshCreation")
-    << " bytes" << std::endl
-    
-    << "MeshStratification " << logger.getNumAllocations("MeshStratification")
-    << " allocations " << logger.getAllocationTotal("MeshStratification")
-    << " bytes" << std::endl
-    
-    << "MeshStratification " << logger.getNumDeallocations("MeshStratification")
-    << " deallocations " << logger.getDeallocationTotal("MeshStratification")
-    << " bytes" << std::endl << std::endl;
-#endif
-
-  assert(0 != _normalizer);
-  const double lengthScale = _normalizer->lengthScale();
-  _normalizer->nondimensionalize(&(*coordinates)[0], coordinates->size(),
-				 lengthScale);
-
-  ALE::SieveBuilder<Mesh>::buildCoordinates(*_mesh, spaceDim, 
-					    &(*coordinates)[0]);
-} // _buildMesh
-
-// ----------------------------------------------------------------------
-// Set vertices and cells for fault mesh.
-void
-pylith::meshio::MeshIO::_buildFaultMesh(const double_array& coordinates,
-					const int numVertices,
-					const int spaceDim,
-					const int_array& cells,
-					const int numCells,
-					const int numCorners,
-					const int firstCell,
-					const int_array& faceCells,
-					const int meshDim,
-					const Obj<Mesh>& fault,
-					Obj<ALE::Mesh>& faultBd)
-{ // _buildFaultMesh
-  MPI_Comm comm = PETSC_COMM_WORLD;
-  int      dim  = meshDim;
-  int      rank;
-
-  assert(!fault.isNull());
-
-  ALE::Obj<sieve_type> sieve = new sieve_type(fault->comm());
-  fault->setDebug(fault->debug());
-  fault->setSieve(sieve);
-
-  MPI_Comm_rank(comm, &rank);
-  // Memory debugging
-  ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
-  logger.setDebug(fault->debug()/2);
-
-  logger.stagePush("FaultCreation");
-  if (!rank) {
-    assert(coordinates.size() == numVertices*spaceDim);
-    assert(cells.size() == numCells*numCorners);
-    {
-      ALE::Obj<ALE::Mesh::sieve_type> s = new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
-
-      ALE::SieveBuilder<ALE::Mesh>::buildTopology(s, meshDim, 
-                                                  numCells, 
-                                                  const_cast<int*>(&cells[0]), 
-                                                  numVertices, 
-                                                  true,
-                                                  numCorners,
-                                                  0,
-                                                  fault->getArrowSection("orientation"),
-                                                  firstCell);
-
-      // Add in cells
-      for(int c = 0; c < numCells; ++c) {
-        s->addArrow(c+firstCell, faceCells[c*2+0]);
-        s->addArrow(c+firstCell, faceCells[c*2+1]);
-      }
-
-      Mesh::renumbering_type& renumbering = fault->getRenumbering();
-      ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering, false);
-      ALE::ISieveConverter::convertOrientation(*s, *sieve, renumbering, fault->getArrowSection("orientation").ptr());
-
-      Obj<ALE::Mesh> tmpMesh = new ALE::Mesh(fault->comm(), dim, fault->debug());
-      faultBd = ALE::Selection<ALE::Mesh>::boundary(tmpMesh);
-    }
-    logger.stagePop();
-    logger.stagePush("FaultStratification");
-    fault->stratify();
-    logger.stagePop();
-  } else {
-    logger.stagePush("FaultStratification");
-    fault->getSieve()->setChart(sieve_type::chart_type());
-    fault->getSieve()->allocate();
-    fault->stratify();
-    logger.stagePop();
-  }
-
-#if defined(ALE_MEM_LOGGING)
-  std::cout
-    << std::endl
-    << "FaultCreation " << logger.getNumAllocations("FaultCreation")
-    << " allocations " << logger.getAllocationTotal("FaultCreation")
-    << " bytes" << std::endl
-    
-    << "FaultCreation " << logger.getNumDeallocations("FaultCreation")
-    << " deallocations " << logger.getDeallocationTotal("FaultCreation")
-    << " bytes" << std::endl
-    
-    << "FaultStratification " << logger.getNumAllocations("FaultStratification")
-    << " allocations " << logger.getAllocationTotal("FaultStratification")
-    << " bytes" << std::endl
-    
-    << "FaultStratification " << logger.getNumDeallocations("FaultStratification")
-    << " deallocations " << logger.getDeallocationTotal("FaultStratification")
-    << " bytes" << std::endl << std::endl;
-#endif
-
-} // _buildFaultMesh
-
-// ----------------------------------------------------------------------
 // Get coordinates of vertices in mesh.
 void
 pylith::meshio::MeshIO::_getVertices(double_array* coordinates,
@@ -353,12 +98,15 @@
   assert(0 != numVertices);
   assert(0 != spaceDim);
   assert(0 != _mesh);
-  assert(!_mesh->isNull());
 
-  const ALE::Obj<Mesh::label_sequence>& vertices = (*_mesh)->depthStratum(0);
+  const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
+  assert(!sieveMesh.isNull());
+
+  const ALE::Obj<SieveMesh::label_sequence>& vertices = 
+    sieveMesh->depthStratum(0);
   assert(!vertices.isNull());
-  const ALE::Obj<real_section_type>& coordsField =
-    (*_mesh)->getRealSection("coordinates");
+  const ALE::Obj<SieveMesh::real_section_type>& coordsField =
+    sieveMesh->getRealSection("coordinates");
   assert(!coordsField.isNull());
 
   *numVertices = vertices->size();
@@ -369,11 +117,13 @@
   const int size = (*numVertices) * (*spaceDim);
   coordinates->resize(size);
 
+  const SieveMesh::label_sequence::iterator verticesEnd = 
+    vertices->end();
   int i = 0;
-  for(Mesh::label_sequence::iterator v_iter = vertices->begin();
-      v_iter != vertices->end();
+  for(SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+      v_iter != verticesEnd;
       ++v_iter) {
-    const real_section_type::value_type *vertexCoords = 
+    const SieveMesh::real_section_type::value_type *vertexCoords = 
       coordsField->restrictPoint(*v_iter);
     for (int iDim=0; iDim < *spaceDim; ++iDim)
       (*coordinates)[i++] = vertexCoords[iDim];
@@ -397,30 +147,34 @@
   assert(0 != numCells);
   assert(0 != meshDim);
   assert(0 != _mesh);
-  assert(!_mesh->isNull());
 
-  const ALE::Obj<sieve_type>& sieve = (*_mesh)->getSieve();
+  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<Mesh::label_sequence>& meshCells = (*_mesh)->heightStratum(0);
+  const ALE::Obj<SieveMesh::label_sequence>& meshCells = 
+    sieveMesh->heightStratum(0);
   assert(!meshCells.isNull());
 
-  *meshDim = (*_mesh)->getDimension();
+  *meshDim = _mesh->dimension();
   *numCells = meshCells->size();
-  *numCorners = (*_mesh)->getNumCellCorners();
+  *numCorners = sieveMesh->getNumCellCorners();
   
-  const ALE::Obj<Mesh::numbering_type>& vNumbering = 
-    (*_mesh)->getFactory()->getLocalNumbering(*_mesh, 0);
+  const ALE::Obj<SieveMesh::numbering_type>& vNumbering = 
+    sieveMesh->getFactory()->getLocalNumbering(sieveMesh, 0);
 
   const int size = (*numCells) * (*numCorners);
   cells->resize(size);
     
-  ALE::ISieveVisitor::PointRetriever<Mesh::sieve_type> pV(sieve->getMaxConeSize());
+  ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> 
+    pV(sieve->getMaxConeSize());
   int i = 0;
-  for(Mesh::label_sequence::iterator e_iter = meshCells->begin();
+  for(SieveMesh::label_sequence::iterator e_iter = meshCells->begin();
       e_iter != meshCells->end();
       ++e_iter) {
     sieve->cone(*e_iter, pV);
-    const Mesh::point_type *cone = pV.getPoints();
+    const SieveMesh::point_type *cone = pV.getPoints();
     for(int p = 0; p < pV.getSize(); ++p, ++i) {
       (*cells)[i] = vNumbering->getIndex(cone[p]);
     }
@@ -434,14 +188,17 @@
 pylith::meshio::MeshIO::_setMaterials(const int_array& materialIds)
 { // _setMaterials
   assert(0 != _mesh);
-  assert(!_mesh->isNull());
 
+  const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
+  assert(!sieveMesh.isNull());
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("MaterialsCreation");
-  const ALE::Obj<Mesh::label_type>& labelMaterials = 
-    (*_mesh)->createLabel("material-id");
-  if (!(*_mesh)->commRank()) {
-    const ALE::Obj<Mesh::label_sequence>& cells = (*_mesh)->heightStratum(0);
+  const ALE::Obj<SieveMesh::label_type>& labelMaterials = 
+    sieveMesh->createLabel("material-id");
+  if (!sieveMesh->commRank()) {
+    const ALE::Obj<SieveMesh::label_sequence>& cells = 
+      sieveMesh->heightStratum(0);
     assert(!cells.isNull());
 
     const int numCells = materialIds.size();
@@ -453,18 +210,22 @@
       throw std::runtime_error(msg.str());
     } // if
     int i = 0;
-    const Mesh::label_sequence::iterator end = cells->end();
+    const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
 
 #ifdef IMESH_NEW_LABELS
-    labelMaterials->setChart((*_mesh)->getSieve()->getChart());
-    for(Mesh::label_sequence::iterator e_iter = cells->begin(); e_iter != end; ++e_iter) {
+    labelMaterials->setChart(sieveMesh->getSieve()->getChart());
+    for(SieveMesh::label_sequence::iterator e_iter = cells->begin();
+	e_iter != cellsEnd;
+	++e_iter) {
       labelMaterials->setConeSize(*e_iter, 1);
     }
     if (cells->size()) {labelMaterials->setSupportSize(0, cells->size());}
     labelMaterials->allocate();
 #endif
-    for(Mesh::label_sequence::iterator e_iter = cells->begin(); e_iter != end; ++e_iter) {
-      (*_mesh)->setValue(labelMaterials, *e_iter, materialIds[i++]);
+    for(SieveMesh::label_sequence::iterator e_iter = cells->begin();
+	e_iter != cellsEnd;
+	++e_iter) {
+      sieveMesh->setValue(labelMaterials, *e_iter, materialIds[i++]);
     }
 #ifdef IMESH_NEW_LABELS
     labelMaterials->recalculateLabel();
@@ -473,8 +234,17 @@
   logger.stagePop();
 
 #if defined(ALE_MEM_LOGGING)
-  std::cout << "MaterialsCreation " << logger.getNumAllocations("MaterialsCreation") << " allocations " << logger.getAllocationTotal("MaterialsCreation") << " bytes" << std::endl;
-  std::cout << "MaterialsCreation " << logger.getNumDeallocations("MaterialsCreation") << " deallocations " << logger.getDeallocationTotal("MaterialsCreation") << " bytes" << std::endl;
+  std::cout 
+    << "MaterialsCreation "
+    << logger.getNumAllocations("MaterialsCreation")
+    << " allocations " << logger.getAllocationTotal("MaterialsCreation")
+    << " bytes"
+    << std::endl
+    << "MaterialsCreation "
+    << logger.getNumDeallocations("MaterialsCreation")
+    << " deallocations " << logger.getDeallocationTotal("MaterialsCreation")
+    << " bytes"
+    << std::endl;
 #endif
 } // _setMaterials
 
@@ -485,25 +255,29 @@
 { // _getMaterials
   assert(0 != materialIds);
   assert(0 != _mesh);
-  assert(!_mesh->isNull());
 
-  const ALE::Obj<Mesh::label_sequence>& cells = (*_mesh)->heightStratum(0);
+  const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
+  assert(!sieveMesh.isNull());
+
+  const ALE::Obj<SieveMesh::label_sequence>& cells = 
+    sieveMesh->heightStratum(0);
   assert(!cells.isNull());
   const int numCells = cells->size();
 
   const int size = numCells;
   materialIds->resize(size);
   
-  const ALE::Obj<Mesh::label_type>& labelMaterials = 
-    (*_mesh)->getLabel("material-id");
+  const ALE::Obj<SieveMesh::label_type>& labelMaterials = 
+    sieveMesh->getLabel("material-id");
   const int idDefault = 0;
-
+  
+  const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
   int i = 0;
-  for(Mesh::label_sequence::iterator e_iter = cells->begin();
-      e_iter != cells->end();
+  for(SieveMesh::label_sequence::iterator e_iter = cells->begin();
+      e_iter != cellsEnd;
       ++e_iter)
     (*materialIds)[i++] = 
-      (*_mesh)->getValue(labelMaterials, *e_iter, idDefault);
+      sieveMesh->getValue(labelMaterials, *e_iter, idDefault);
 } // _getMaterials
 
 // ----------------------------------------------------------------------
@@ -513,32 +287,44 @@
 				  const GroupPtType type,
 				  const int_array& points)
 { // _setGroup
+  typedef SieveMesh::int_section_type::chart_type chart_type;
+
   assert(0 != _mesh);
-  assert(!_mesh->isNull());
 
+  const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
+  assert(!sieveMesh.isNull());
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("GroupCreation");
-  const ALE::Obj<int_section_type>& groupField = (*_mesh)->getIntSection(name);
+  const ALE::Obj<SieveMesh::int_section_type>& groupField = 
+    sieveMesh->getIntSection(name);
   assert(!groupField.isNull());
 
   const int numPoints   = points.size();
-  const int numVertices = (*_mesh)->depthStratum(0)->size();
-  const int numCells    = (*_mesh)->heightStratum(0)->size();
+  const int numVertices = sieveMesh->depthStratum(0)->size();
+  const int numCells    = sieveMesh->heightStratum(0)->size();
   if (CELL == type) {
-    groupField->setChart(int_section_type::chart_type(0,numCells));
+    groupField->setChart(chart_type(0,numCells));
     for(int i=0; i < numPoints; ++i)
       groupField->setFiberDimension(points[i], 1);
   } else if (VERTEX == type) {
-    groupField->setChart(int_section_type::chart_type(numCells,numCells+numVertices));
+    groupField->setChart(chart_type(numCells, numCells+numVertices));
     for(int i=0; i < numPoints; ++i)
       groupField->setFiberDimension(numCells+points[i], 1);
   } // if/else
-  (*_mesh)->allocate(groupField);
+  sieveMesh->allocate(groupField);
   logger.stagePop();
 
 #if defined(ALE_MEM_LOGGING)
-  std::cout << "GroupCreation " << logger.getNumAllocations("GroupCreation") << " allocations " << logger.getAllocationTotal("GroupCreation") << " bytes" << std::endl;
-  std::cout << "GroupCreation " << logger.getNumDeallocations("GroupCreation") << " deallocations " << logger.getDeallocationTotal("GroupCreation") << " bytes" << std::endl;
+  std::cout
+    << "GroupCreation " << logger.getNumAllocations("GroupCreation")
+    << " allocations " << logger.getAllocationTotal("GroupCreation")
+    << " bytes"
+    << std::endl
+    << "GroupCreation " << logger.getNumDeallocations("GroupCreation")
+    << " deallocations " << logger.getDeallocationTotal("GroupCreation")
+    << " bytes"
+    << std::endl;
 #endif
 } // _setGroup
 
@@ -548,36 +334,42 @@
 pylith::meshio::MeshIO::_distributeGroups()
 { // _distributeGroups
   assert(0 != _mesh);
-  assert(!_mesh->isNull());
 
-  if (!(*_mesh)->commRank()) {
+  const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
+  assert(!sieveMesh.isNull());
+
+  if (!sieveMesh->commRank()) {
     const ALE::Obj<std::set<std::string> >& sectionNames = 
-      (*_mesh)->getIntSections();
+      sieveMesh->getIntSections();
     int numGroups = sectionNames->size();
 
-    MPI_Bcast(&numGroups, 1, MPI_INT, 0, (*_mesh)->comm());
+    MPI_Bcast(&numGroups, 1, MPI_INT, 0, sieveMesh->comm());
+
+    const std::set<std::string>::const_iterator namesEnd = sectionNames->end();
     for (std::set<std::string>::const_iterator name=sectionNames->begin();
-         name != sectionNames->end(); ++name) {
+         name != namesEnd;
+	 ++name) {
       int len = name->size();
       
-      MPI_Bcast(&len, 1, MPI_INT, 0, (*_mesh)->comm());
-      MPI_Bcast((void *) name->c_str(), len, MPI_CHAR, 0, (*_mesh)->comm());
+      MPI_Bcast(&len, 1, MPI_INT, 0, sieveMesh->comm());
+      MPI_Bcast((void *) name->c_str(), len, MPI_CHAR, 0, sieveMesh->comm());
     }
   } else {
     int numGroups;
 
-    MPI_Bcast(&numGroups, 1, MPI_INT, 0, (*_mesh)->comm());
+    MPI_Bcast(&numGroups, 1, MPI_INT, 0, sieveMesh->comm());
     for(int g = 0; g < numGroups; ++g) {
       char *name;
       int   len;
 
-      MPI_Bcast(&len, 1, MPI_INT, 0, (*_mesh)->comm());
+      MPI_Bcast(&len, 1, MPI_INT, 0, sieveMesh->comm());
       name = new char[len+1];
-      MPI_Bcast(name, len, MPI_CHAR, 0, (*_mesh)->comm());
+      MPI_Bcast(name, len, MPI_CHAR, 0, sieveMesh->comm());
       name[len] = 0;
-      const ALE::Obj<int_section_type>& groupField = (*_mesh)->getIntSection(name);
+      const ALE::Obj<SieveMesh::int_section_type>& groupField = 
+	sieveMesh->getIntSection(name);
       assert(!groupField.isNull());
-      (*_mesh)->allocate(groupField);
+      sieveMesh->allocate(groupField);
       delete [] name;
     } // for
   } // if/else
@@ -590,16 +382,19 @@
 { // _getGroups
   assert(0 != names);
   assert(0 != _mesh);
-  assert(!_mesh->isNull());
-  
+
+  const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
+  assert(!sieveMesh.isNull());
+
   const ALE::Obj<std::set<std::string> >& sectionNames = 
-    (*_mesh)->getIntSections();
+    sieveMesh->getIntSections();
   
   const int numGroups = sectionNames->size();
   names->resize(numGroups);
+  const std::set<std::string>::const_iterator namesEnd = sectionNames->end();
   int i=0;
   for (std::set<std::string>::const_iterator name=sectionNames->begin();
-       name != sectionNames->end();
+       name != namesEnd;
        ++name)
     (*names)[i++] = *name;
 } // _getGroups
@@ -611,39 +406,52 @@
 				  GroupPtType* type,
 				  const char *name) const
 { // _getGroup
+  typedef SieveMesh::int_section_type::chart_type chart_type;
+
   assert(0 != points);
   assert(0 != type);
   assert(0 != _mesh);
-  assert(!_mesh->isNull());
 
-  const ALE::Obj<int_section_type>& groupField = (*_mesh)->getIntSection(name);
+  const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
+  assert(!sieveMesh.isNull());
+
+  const ALE::Obj<SieveMesh::int_section_type>& groupField = 
+    sieveMesh->getIntSection(name);
   assert(!groupField.isNull());
-  const int_section_type::chart_type& chart = groupField->getChart();
-  Mesh::point_type firstPoint;
-  for(int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
-    if (groupField->getFiberDimension(*c_iter)) {firstPoint = *c_iter; break;}
+  const chart_type& chart = groupField->getChart();
+  SieveMesh::point_type firstPoint;
+  chart_type::const_iterator chartEnd = chart.end();
+  for(chart_type::const_iterator c_iter=chart.begin();
+      c_iter != chartEnd;
+      ++c_iter) {
+    if (groupField->getFiberDimension(*c_iter)) {
+      firstPoint = *c_iter;
+      break;
+    } // if
   }
-  ALE::Obj<Mesh::numbering_type> numbering;
+  ALE::Obj<SieveMesh::numbering_type> numbering;
 
-  if ((*_mesh)->height(firstPoint) == 0) {
+  if (sieveMesh->height(firstPoint) == 0) {
     *type = CELL;
-    numbering = (*_mesh)->getFactory()->getNumbering(*_mesh, 
-                                                    (*_mesh)->depth());
+    numbering = sieveMesh->getFactory()->getNumbering(sieveMesh, 
+						      sieveMesh->depth());
   } else {
     *type = VERTEX;
-    numbering = (*_mesh)->getFactory()->getNumbering(*_mesh, 0);
+    numbering = sieveMesh->getFactory()->getNumbering(sieveMesh, 0);
   } // if/else
   const int numPoints = groupField->size();
   points->resize(numPoints);
   int i = 0;
 
-  for(int_section_type::chart_type::const_iterator c_iter = chart.begin();
-      c_iter != chart.end();
+  for(chart_type::const_iterator c_iter=chart.begin();
+      c_iter != chartEnd;
       ++c_iter) {
     assert(!numbering.isNull());
-    if (groupField->getFiberDimension(*c_iter)) (*points)[i++] = numbering->getIndex(*c_iter);
+    if (groupField->getFiberDimension(*c_iter)) (*points)[i++] = 
+      numbering->getIndex(*c_iter);
   } // for
 } // _getGroup
 
 
-// End of file 
+// End of file
+

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIO.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIO.hh	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIO.hh	2009-01-24 00:22:09 UTC (rev 13936)
@@ -13,24 +13,28 @@
 #if !defined(pylith_meshio_meshio_hh)
 #define pylith_meshio_meshio_hh
 
+// Include directives ---------------------------------------------------
 #include "pylith/utils/arrayfwd.hh" // USES double_array, int_array,
                                     // string_vector
 
-#include "pylith/utils/sievetypes.hh" // USES Obj, PETSc Mesh
-#include "Mesh.hh" // USES ALE::Mesh
-
+// Forward declarations -------------------------------------------------
 namespace pylith {
   namespace meshio {
     class MeshIO;
   } // meshio
+
+  namespace topology {
+    class Mesh; // USES Mesh
+  } // topology
 } // pylith
 
 namespace spatialdata {
   namespace units {
-    class Nondimensional;
+    class Nondimensional; // USES Nondimensional
   } // units
 } // spatialdata
 
+// MeshIO ---------------------------------------------------------------
 class pylith::meshio::MeshIO
 { // MeshIO
 
@@ -38,7 +42,10 @@
 public :
 
   /// Type of points in a group.
-  typedef enum { VERTEX, CELL } GroupPtType;
+  enum GroupPtType {
+    VERTEX=0,
+    CELL=1,
+  }; // GroupPtType
 
 // PUBLIC MEMBERS ///////////////////////////////////////////////////////
 public :
@@ -46,7 +53,8 @@
   MeshIO(void);
 
   /// Destructor
-  virtual ~MeshIO(void);
+  virtual
+  ~MeshIO(void);
 
   /** Set debug flag for mesh.
    *
@@ -82,15 +90,15 @@
 
   /** Read mesh from file.
    *
-   * @param mesh Pointer to PETSc mesh object
+   * @param mesh PyLith finite-element mesh.
    */
-  void read(ALE::Obj<Mesh>* mesh);
+  void read(topology::Mesh* mesh);
 
   /** Write mesh to file.
    *
-   * @param mesh Pointer to PETSc mesh object.
+   * @param mesh PyLith finite-element mesh.
    */
-  void write(ALE::Obj<Mesh>* mesh);
+  void write(topology::Mesh* const mesh);
 
 // PROTECTED MEMBERS ////////////////////////////////////////////////////
 protected :
@@ -109,40 +117,6 @@
    */
   int getMeshDim(void) const;
 
-  /** Build mesh topology and set vertex coordinates.
-   *
-   * All mesh information must use zero based indices. In other words,
-   * the lowest index MUST be 0 not 1.
-   *
-   * @param coordinates Array of coordinates of vertices
-   * @param numVertices Number of vertices
-   * @param spaceDim Dimension of vector space for vertex coordinates
-   * @param cells Array of indices of vertices in cells (first index is 0 or 1)
-   * @param numCells Number of cells
-   * @param numCorners Number of vertices per cell
-   * @param meshDim Dimension of cells in mesh
-   */
-  void _buildMesh(double_array* coordinates,
-		  const int numVertices,
-		  const int spaceDim,
-		  const int_array& cells,
-		  const int numCells,
-		  const int numCorners,
-		  const int meshDim);
-
-  static void
-  _buildFaultMesh(const double_array& coordinates,
-                  const int numVertices,
-                  const int spaceDim,
-                  const int_array& cells,
-                  const int numCells,
-                  const int numCorners,
-                  const int firstCell,
-                  const int_array& faceCells,
-                  const int meshDim,
-                  const Obj<Mesh>& fault,
-                  Obj<ALE::Mesh>& faultBd);
-
   /** Get information about vertices in mesh.
    *
    * Method caller is responsible for memory management.
@@ -218,18 +192,17 @@
 		 GroupPtType* type,
 		 const char *name) const;
 
-  /** Create empty groups on other processes
-   */
+  /// Create empty groups on other processes
   void _distributeGroups();
 
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
 
-  ALE::Obj<Mesh>* _mesh; ///< Pointer to PETSc mesh object
-  spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer
+  topology::Mesh* _mesh; ///< Pointer to finite-element mesh.
+  spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer.
 
-  bool _debug; ///< True to turn of mesh debugging output
-  bool _interpolate; ///< True if building intermediate topology elements
+  bool _debug; ///< True to turn of mesh debugging output.
+  bool _interpolate; ///< True if building intermediate topology elements.
 
 }; // MeshIO
 
@@ -237,4 +210,5 @@
 
 #endif // pylith_meshio_meshio_hh
 
+
 // End of file 

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOAscii.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOAscii.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOAscii.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -14,13 +14,14 @@
 
 #include "MeshIOAscii.hh" // implementation of class methods
 
+#include "MeshBuilder.hh" // USES MeshBuilder
+#include "pylith/topology/Mesh.hh" // USES Mesh
+
 #include "pylith/utils/array.hh" // USES double_array, int_array, string_vector
 #include "spatialdata/utils/LineParser.hh" // USES LineParser
 
 #include "journal/info.h" // USES journal::info_t
 
-#include <petsc.h> // USES MPI
-
 #include <iomanip> // USES setw(), setiosflags(), resetiosflags()
 #include <strings.h> // USES strcasecmp()
 #include <cassert> // USES assert()
@@ -29,8 +30,10 @@
 #include <sstream> // USES std::ostringstream
 
 // ----------------------------------------------------------------------
-const char* pylith::meshio::MeshIOAscii::groupTypeNames[] = 
-  {"vertices", "cells"};
+const char* pylith::meshio::MeshIOAscii::groupTypeNames[] = {
+  "vertices",
+  "cells",
+};
 
 // ----------------------------------------------------------------------
 // Constructor
@@ -51,8 +54,10 @@
 void
 pylith::meshio::MeshIOAscii::_read(void)
 { // _read
-  MPI_Comm comm = PETSC_COMM_WORLD;
-  int rank;
+  assert(0 != _normalizer);
+
+  MPI_Comm comm = _mesh->comm();
+  int rank = 0;
   int meshDim = 0;
   int spaceDim = 0;
   int numVertices = 0;
@@ -63,7 +68,7 @@
   int_array materialIds;
 
   MPI_Comm_rank(comm, &rank);
-  if (!rank) {
+  if (0 == rank) {
     std::ifstream filein(_filename.c_str());
     if (!filein.is_open() || !filein.good()) {
       std::ostringstream msg;
@@ -74,11 +79,11 @@
 
     spatialdata::utils::LineParser parser(filein, "//");
     parser.eatwhitespace(true);
-
+    
     std::string token;
     std::istringstream buffer;
     const int maxIgnore = 1024;
-
+    
     buffer.str(parser.next());
     buffer >> token;
     if (0 != strcasecmp(token.c_str(), "mesh")) {
@@ -134,8 +139,9 @@
 
 	if (readDim && readCells && readVertices && !builtMesh) {
 	  // Can now build mesh
-	  _buildMesh(&coordinates, numVertices, spaceDim,
-		     cells, numCells, numCorners, meshDim);
+	  MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
+				 cells, numCells, numCorners, meshDim,
+				 _interpolate, *_normalizer);
 	  _setMaterials(materialIds);
 	  builtMesh = true;
 	} // if
@@ -161,8 +167,9 @@
     } // catch
     filein.close();
   } else {
-    _buildMesh(&coordinates, numVertices, spaceDim,
-               cells, numCells, numCorners, meshDim);
+    MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
+			   cells, numCells, numCorners, meshDim,
+			   _interpolate, *_normalizer);
     _setMaterials(materialIds);
   } // if/else
   _distributeGroups();

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOAscii.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOAscii.hh	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOAscii.hh	2009-01-24 00:22:09 UTC (rev 13936)
@@ -10,17 +10,27 @@
 // ======================================================================
 //
 
+/**
+ * @file libsrc/meshio/MeshIOAscii.hh
+ *
+ * @brief C++ input/output manager for PyLith ASCII mesh files.
+ */
+
 #if !defined(pylith_meshio_meshioascii_hh)
 #define pylith_meshio_meshioascii_hh
 
+// Include directives ---------------------------------------------------
+#include "MeshIO.hh" // ISA MeshIO
+
 #include <iosfwd> // USES std::istream, std::ostream
 #include <string> // HASA std::string
 
-#include "MeshIO.hh"
-
+// Forward declarations -------------------------------------------------
 namespace pylith {
   namespace meshio {
     class MeshIOAscii;
+
+    class TestMeshIOAscii; // unit testing
   } // meshio
 } // pylith
 
@@ -30,12 +40,13 @@
   } // utils
 } // pylith
 
+// MeshIOAscii ----------------------------------------------------------
 class pylith::meshio::MeshIOAscii : public MeshIO
 { // MeshIOAscii
+  friend class TestMeshIOAscii; // unit testing
 
 // PUBLIC METHODS ///////////////////////////////////////////////////////
 public :
-  static const char *groupTypeNames[];
 
   /// Constructor
   MeshIOAscii(void);
@@ -132,6 +143,9 @@
   std::string _filename; ///< Name of file
   bool _useIndexZero; ///< Flag indicating if indicates start at 0 (T) or 1 (F)
 
+  static
+  const char *groupTypeNames[]; ///< Types of mesh groups.
+
 }; // MeshIOAscii
 
 #include "MeshIOAscii.icc" // inline methods

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOCubit.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOCubit.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOCubit.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -14,6 +14,8 @@
 
 #include "MeshIOCubit.hh" // implementation of class methods
 
+#include "MeshBuilder.hh" // USES MeshBuilder
+
 #include "pylith/utils/array.hh" // USES double_array, int_array, string_vector
 
 #include "journal/info.h" // USES journal::info_t
@@ -45,7 +47,7 @@
 pylith::meshio::MeshIOCubit::_read(void)
 { // _read
   MPI_Comm comm = PETSC_COMM_WORLD;
-  int rank;
+  int rank = 0;
   int meshDim = 0;
   int spaceDim = 0;
   int numVertices = 0;
@@ -56,7 +58,7 @@
   int_array materialIds;
 
   MPI_Comm_rank(comm, &rank);
-  if (!rank) {
+  if (0 == rank) {
     try {
       NcFile ncfile(_filename.c_str());
       if (!ncfile.is_valid()) {
@@ -75,10 +77,11 @@
       _readVertices(ncfile, &coordinates, &numVertices, &spaceDim);
       _readCells(ncfile, &cells, &materialIds, &numCells, &numCorners);
       _orientCells(&cells, numCells, numCorners, meshDim);
-      _buildMesh(&coordinates, numVertices, spaceDim,
-                 cells, numCells, numCorners, meshDim);
+      MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
+			     cells, numCells, numCorners, meshDim,
+			     _interpolate, *_normalizer);
       _setMaterials(materialIds);
-
+      
       _readGroups(ncfile);
     } catch (std::exception& err) {
       std::ostringstream msg;
@@ -92,8 +95,9 @@
       throw std::runtime_error(msg.str());
     } // try/catch
   } else {
-    _buildMesh(&coordinates, numVertices, spaceDim,
-               cells, numCells, numCorners, meshDim);
+    MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
+			   cells, numCells, numCorners, meshDim,
+			   _interpolate, *_normalizer);
     _setMaterials(materialIds);
   }
   _distributeGroups();
@@ -219,7 +223,8 @@
       std::ostringstream msg;
       msg << "All materials must have the same number of vertices per cell.\n"
 	  << "Expected " << *numCorners << " vertices per cell, but block "
-	  << blockIds[iMaterial] << " has " << num_nod_per_el->size() << " vertices.";
+	  << blockIds[iMaterial] << " has " << num_nod_per_el->size()
+	  << " vertices.";
       throw std::runtime_error(msg.str());
     } // if
 

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOCubit.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOCubit.hh	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOCubit.hh	2009-01-24 00:22:09 UTC (rev 13936)
@@ -10,13 +10,21 @@
 // ======================================================================
 //
 
+/**
+ * @file pylith/meshio/MeshIOCubit.hh
+ *
+ * @brief C++ input/output manager for CUBIT Exodus II files.
+ */
+
 #if !defined(pylith_meshio_meshiocubit_hh)
 #define pylith_meshio_meshiocubit_hh
 
+// Include directives ---------------------------------------------------
+#include "MeshIO.hh" // ISA MeshIO
+
 #include <string> // HASA std::string
 
-#include "MeshIO.hh"
-
+// Forward declarations -------------------------------------------------
 namespace pylith {
   namespace meshio {
     class MeshIOCubit;
@@ -27,9 +35,10 @@
 
 class NcFile; // netcdf file
 
+// MeshIOCubit ----------------------------------------------------------
 class pylith::meshio::MeshIOCubit : public MeshIO
 { // MeshIOCubit
-  friend class TestMeshIOCubit;
+  friend class TestMeshIOCubit; // unit testing
 
 // PUBLIC METHODS ///////////////////////////////////////////////////////
 public :

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOLagrit.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOLagrit.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOLagrit.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -18,13 +18,13 @@
 #include "GMVFileBinary.hh" // USES GMVFileBinary
 #include "PsetFileAscii.hh" // USES PsetFileAscii
 #include "PsetFileBinary.hh" // USES PsetFileBinary
+#include "MeshBuilder.hh" // USES MeshBuilder
 
+#include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/utils/array.hh" // USES double_array, int_array
 
-#include <petsc.h> // USES MPI_Comm
-
 #include <cassert> // USES assert()
-#include <stdexcept> // TEMPORARY
+#include <stdexcept> // USES std::runtime_error()
 
 // ----------------------------------------------------------------------
 // Constructor
@@ -43,151 +43,12 @@
 { // destructor
 } // destructor
 
-void
-pylith::meshio::MeshIOLagrit::readFault(const std::string filename, const Obj<Mesh>& mesh, const Obj<Mesh>& fault, Obj<ALE::Mesh>& faultBd) {
-  int faultDim = 2;
-  int fSpaceDim = 0;
-  int numFVertices = 0;
-  int numFCells = 0;
-  int numFCorners = 0;
-  double_array fCoordinates;
-  int_array fCells;
-  int_array fMaterialIds;
-  int_array faceCells;
-  int_array vertexIDs;
-
-  if (!fault->commRank()) {
-    std::ifstream fin(filename.c_str(), std::ios::in);
-    if (!(fin.is_open() && fin.good())) {
-      std::ostringstream msg;
-      msg << "Could not open ASCII INP file '" << filename << "' for reading.";
-      throw std::runtime_error(msg.str());
-    } // if
-    int numNodeData, numCellData, numModelData;
-
-    fSpaceDim = 3;
-    // Section 1: <num_nodes> <num_cells> <num_ndata> <num_cdata> <num_mdata>
-    fin >> numFVertices;
-    fin >> numFCells;
-    fin >> numNodeData;
-    fin >> numCellData;
-    fin >> numModelData;
-    // Section 2: <node_id 1> <x> <y> <z>
-    fCoordinates.resize(numFVertices * fSpaceDim);
-    for(int v = 0; v < numFVertices; ++v) {
-      int id;
-
-      fin >> id;
-      for(int d = 0; d < fSpaceDim; ++d) {
-        fin >> fCoordinates[v*fSpaceDim + d];
-      }
-    }
-    // Section 3: <cell_id 1> <mat_id> <cell_type> <cell_vert 1> ... <cell_vert n> 
-    fMaterialIds.resize(numFCells);
-    for(int c = 0; c < numFCells; ++c) {
-      std::string cellType;
-      int         cellID;
-
-      fin >> cellID;
-      fin >> fMaterialIds[c];
-      fin >> cellType;
-      if (cellType == "tri") {
-        numFCorners = 3;
-      } else if (cellType == "quad") {
-        numFCorners = 4;
-      } else {
-        std::ostringstream msg;
-        msg << "Unknown cell type " << cellType << "while reading INP file.";
-        throw std::runtime_error(msg.str());
-      }
-      if (c == 0) fCells.resize(numFCells*numFCorners);
-      for(int i = 0; i < numFCorners; ++i) fin >> fCells[c*numFCorners+i];
-    }
-    // Section 4: <num_comp for node data> <size comp 1> <size comp 2>...<size comp n>
-    int numComponents, totalSize;
-
-    fin >> numComponents;
-    totalSize = 0;
-    for(int i = 0; i < numComponents; ++i) {
-      int compSize;
-      fin >> compSize;
-      totalSize += compSize;
-    }
-    // Section 5: <node_comp_label 1> , <units_label 1>
-    for(int c = 0; c < numComponents; ++c) {
-      std::string label, typeName;
-
-      fin >> label;
-      fin >> typeName;
-    }
-    // Section 6: <node_id 1> <node_data 1> ... <node_data num_ndata>
-    vertexIDs.resize(numFVertices);
-    for(int v = 0; v < numFVertices; ++v) {
-      int id;
-      double dummy;
-
-      fin >> id;
-      fin >> vertexIDs[v]; // global node number
-      fin >> dummy; // components of normal at vertices
-      fin >> dummy;
-      fin >> dummy;
-    }
-    // Section 7: <num_comp for cell's data> <size comp 1> <size comp 2>...<size comp n>
-    fin >> numComponents;
-    totalSize = 0;
-    for(int i = 0; i < numComponents; ++i) {
-      int compSize;
-      fin >> compSize;
-      totalSize += compSize;
-    }
-    // Section 8: <cell_component_label 1> , <units_label 1> 
-    for(int c = 0; c < numComponents; ++c) {
-      std::string label, typeName;
-
-      fin >> label;
-      fin >> typeName;
-    }
-    // Section 9: <cell_id 1> <cell_data 1> ... <cell_data num_cdata> 
-    faceCells.resize(numFCells*2);
-    for(int c = 0; c < numFCells; ++c) {
-      int id, faultId;
-
-      fin >> id;
-      fin >> faceCells[c*2+0]; // Cell numbers in global mesh on either side of fault
-      fin >> faceCells[c*2+1];
-      fin >> faultId;
-    }
-
-    // Determine the number of cells
-    //   Only do this once since we add cohesive cells after that
-    static int numCells = -1;
-
-    if (numCells == -1) {numCells = mesh->heightStratum(0)->size();}
-
-    // Renumber vertices and use zero based indices
-    // UCD file has one-based indices for both vertexIDs and fCells
-    //   Also, vertex numbers are offset by the number of cells
-    for(int c = 0; c < numFCells; ++c)
-      for(int corner = 0; corner < numFCorners; ++corner)
-        fCells[c*numFCorners+corner] = 
-	  vertexIDs[fCells[c*numFCorners+corner]-1] - 1 + numCells;
-
-    // Switch to zero based index for global cell numbering
-    for (int c=0; c < numFCells; ++c)
-      for (int i=0; i < 2; ++i)
-	faceCells[c*2+i] -= 1;
-  }
-
-  const int firstFaultCell = mesh->getSieve()->getBaseSize() + mesh->getSieve()->getCapSize();
-  _buildFaultMesh(fCoordinates, numFVertices, fSpaceDim, fCells, numFCells, numFCorners, firstFaultCell, faceCells, faultDim, fault, faultBd);
-}
-
 // ----------------------------------------------------------------------
 // Unpickle mesh
 void
 pylith::meshio::MeshIOLagrit::_read(void)
 { // _read
-  MPI_Comm comm = PETSC_COMM_WORLD;
+  MPI_Comm comm = _mesh->comm();
   int rank;
   int meshDim = 0;
   int spaceDim = 0;
@@ -212,11 +73,12 @@
       _orientCellsBinary(&cells, numCells, numCorners, meshDim);
     } // if/else
   }
-  _buildMesh(&coordinates, numVertices, spaceDim,
-             cells, numCells, numCorners, meshDim);
+  MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
+			 cells, numCells, numCorners, meshDim,
+			 _interpolate, *_normalizer);
   _setMaterials(materialIds);
 
-  if (!rank) {
+  if (0 == rank) {
     std::vector<PsetFile::Pset> groups;
     if (PsetFile::isAscii(_filenamePset.c_str())) {
       PsetFileAscii filein(_filenamePset.c_str());

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOLagrit.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOLagrit.hh	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/MeshIOLagrit.hh	2009-01-24 00:22:09 UTC (rev 13936)
@@ -10,13 +10,21 @@
 // ======================================================================
 //
 
+/**
+ * @file pylith/meshio/MeshIOLagrit.hh
+ *
+ * @brief Input/output manager for LaGriT GMV and Pset files.
+ */
+
 #if !defined(pylith_meshio_meshiolagrit_hh)
 #define pylith_meshio_meshiolagrit_hh
 
+// Include directives ---------------------------------------------------
 #include "MeshIO.hh" // ISA MeshIO
 
 #include <string> // HASA std::string
 
+// Forward declarations -------------------------------------------------
 namespace pylith {
   namespace meshio {
     class MeshIOLagrit;
@@ -25,6 +33,7 @@
   } // meshio
 } // pylith
 
+// MeshIOLagrit ---------------------------------------------------------
 class pylith::meshio::MeshIOLagrit : public MeshIO
 { // MeshIOLagrit
   friend class TestMeshIOLagrit; // unit testing
@@ -111,13 +120,6 @@
    */
   bool isRecordHeader32Bit(void) const;
 
-  // TEMPORARY
-  static
-  void readFault(const std::string filename,
-		 const Obj<Mesh>& mesh,
-		 const Obj<Mesh>& fault,
-		 Obj<ALE::Mesh>& faultBd);
-
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :
 

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/Makefile.am	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/Makefile.am	2009-01-24 00:22:09 UTC (rev 13936)
@@ -21,38 +21,40 @@
 
 # Primary source files
 testmeshio_SOURCES = \
-	TestCellFilterAvg.cc \
-	TestDataWriterVTK.cc \
-	TestDataWriterVTKMesh.cc \
-	TestDataWriterVTKMeshLine2.cc \
-	TestDataWriterVTKMeshTri3.cc \
-	TestDataWriterVTKMeshQuad4.cc \
-	TestDataWriterVTKMeshTet4.cc \
-	TestDataWriterVTKMeshHex8.cc \
-	TestDataWriterVTKSubMeshLine2.cc \
-	TestDataWriterVTKSubMeshTri3.cc \
-	TestDataWriterVTKSubMeshQuad4.cc \
-	TestDataWriterVTKSubMeshTet4.cc \
-	TestDataWriterVTKSubMeshHex8.cc \
-	TestDataWriterVTKFaultMesh.cc \
-	TestDataWriterVTKFaultMeshTri3.cc \
-	TestDataWriterVTKFaultMeshQuad4.cc \
-	TestDataWriterVTKFaultMeshTet4.cc \
-	TestDataWriterVTKFaultMeshHex8.cc \
-	TestDataWriterVTKBCMesh.cc \
-	TestDataWriterVTKBCMeshTri3.cc \
-	TestDataWriterVTKBCMeshQuad4.cc \
-	TestDataWriterVTKBCMeshTet4.cc \
-	TestDataWriterVTKBCMeshHex8.cc \
 	TestMeshIO.cc \
 	TestMeshIOAscii.cc \
 	TestMeshIOLagrit.cc \
-	TestOutputManager.cc \
-	TestOutputSolnSubset.cc \
-	TestVertexFilterVecNorm.cc \
 	test_meshio.cc
 
+#	TestCellFilterAvg.cc \
+# 	TestDataWriterVTK.cc \
+# 	TestDataWriterVTKMesh.cc \
+# 	TestDataWriterVTKMeshLine2.cc \
+# 	TestDataWriterVTKMeshTri3.cc \
+# 	TestDataWriterVTKMeshQuad4.cc \
+# 	TestDataWriterVTKMeshTet4.cc \
+# 	TestDataWriterVTKMeshHex8.cc \
+# 	TestDataWriterVTKSubMeshLine2.cc \
+# 	TestDataWriterVTKSubMeshTri3.cc \
+# 	TestDataWriterVTKSubMeshQuad4.cc \
+# 	TestDataWriterVTKSubMeshTet4.cc \
+# 	TestDataWriterVTKSubMeshHex8.cc \
+# 	TestDataWriterVTKFaultMesh.cc \
+# 	TestDataWriterVTKFaultMeshTri3.cc \
+# 	TestDataWriterVTKFaultMeshQuad4.cc \
+# 	TestDataWriterVTKFaultMeshTet4.cc \
+# 	TestDataWriterVTKFaultMeshHex8.cc \
+# 	TestDataWriterVTKBCMesh.cc \
+# 	TestDataWriterVTKBCMeshTri3.cc \
+# 	TestDataWriterVTKBCMeshQuad4.cc \
+# 	TestDataWriterVTKBCMeshTet4.cc \
+# 	TestDataWriterVTKBCMeshHex8.cc
 
+# 	TestOutputManager.cc \
+# 	TestOutputSolnSubset.cc \
+# 	TestVertexFilterVecNorm.cc
+
+
 noinst_HEADERS = \
 	TestCellFilterAvg.hh \
 	TestDataWriterVTK.hh \

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestCellFilterAvg.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestCellFilterAvg.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestCellFilterAvg.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -17,6 +17,7 @@
 #include "pylith/meshio/CellFilterAvg.hh"
 
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+#include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/feassemble/Quadrature2D.hh" // USES Quadrature
 
 // ----------------------------------------------------------------------
@@ -76,25 +77,29 @@
   };
 
   MeshIOAscii iohandler;
-  ALE::Obj<Mesh> mesh;
+  topology::Mesh(PETSC_COMM_WORLD, cellDim);
   iohandler.filename(filename);
   iohandler.read(&mesh);
   CPPUNIT_ASSERT(!mesh.isNull());
 
+  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+  CPPUNIT_ASSERT(!sieveMesh.isNull());
+
   // Set cell field
-  const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
-  const Mesh::label_sequence::iterator cellsEnd = cells->end();
+  const ALE::Obj<SieveMesh::label_sequence>& cells = 
+    sieveMesh->heightStratum(0);
+  const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
 
-  ALE::Obj<real_section_type> field = 
-    new real_section_type(mesh->comm(), mesh->debug());
-  field->setChart(real_section_type::chart_type(0, cells->size()));
+  ALE::Obj<SieveMesh::real_section_type> field = 
+    new SieveMesh::real_section_type(mesh.comm(), mesh.debug());
+  field->setChart(SieveMesh::real_section_type::chart_type(0, cells->size()));
   field->setFiberDimension(cells, fiberDim);
-  mesh->allocate(field);
+  sieveMesh->allocate(field);
 
   CPPUNIT_ASSERT_EQUAL(ncells, int(cells->size()));
 
   int ipt = 0;
-  for (Mesh::label_sequence::iterator c_iter=cells->begin();
+  for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
        c_iter != cellsEnd;
        ++c_iter, ++ipt) {
     const double* values = &fieldValues[ipt*fiberDim];
@@ -110,11 +115,11 @@
 
   VectorFieldEnum fieldTypeF = SCALAR_FIELD;
   const ALE::Obj<real_section_type>& fieldF =
-    filter.filter(&fieldTypeF, field, mesh);
+    filter.filter(&fieldTypeF, field, meshMesh);
 
   CPPUNIT_ASSERT_EQUAL(fieldTypeE, fieldTypeF);
   ipt = 0;
-  for (Mesh::label_sequence::iterator c_iter=cells->begin();
+  for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
        c_iter != cellsEnd;
        ++c_iter, ++ipt) {
     CPPUNIT_ASSERT_EQUAL(fiberDimE, fieldF->getFiberDimension(*c_iter));

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestCellFilterAvg.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestCellFilterAvg.hh	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestCellFilterAvg.hh	2009-01-24 00:22:09 UTC (rev 13936)
@@ -21,16 +21,17 @@
 #if !defined(pylith_meshio_testcellfilteravg_hh)
 #define pylith_meshio_testcellfilteravg_hh
 
+// Include directives ---------------------------------------------------
 #include <cppunit/extensions/HelperMacros.h>
 
-/// Namespace for pylith package
+// Forward declarations -------------------------------------------------
 namespace pylith {
   namespace meshio {
     class TestCellFilterAvg;
   } // meshio
 } // pylith
 
-/// C++ unit testing for CellFilterAvg
+// TestCellFilterAvg ----------------------------------------------------
 class pylith::meshio::TestCellFilterAvg : public CppUnit::TestFixture
 { // class TestCellFilterAvg
 
@@ -55,4 +56,5 @@
 
 #endif // pylith_meshio_testcellfilteravg_hh
 
+
 // End of file 

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIO.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIO.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIO.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -14,8 +14,9 @@
 
 #include "TestMeshIO.hh" // Implementation of class methods
 
+#include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/meshio/MeshIO.hh" // USES MeshIO
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+
 #include "pylith/utils/array.hh" // USES int_array
 
 #include "data/MeshData.hh"
@@ -25,7 +26,7 @@
 
 // ----------------------------------------------------------------------
 // Get simple mesh for testing I/O.
-ALE::Obj<pylith::Mesh>*
+pylith::topology::Mesh*
 pylith::meshio::TestMeshIO::_createMesh(const MeshData& data)
 { // _createMesh
   // buildTopology() requires zero based index
@@ -41,45 +42,49 @@
     CPPUNIT_ASSERT(0 != data.groupTypes);
   } // if
 
-  ALE::Obj<Mesh>* mesh = new ALE::Obj<Mesh>;
+  topology::Mesh* mesh = new topology::Mesh(PETSC_COMM_WORLD, data.cellDim);
   CPPUNIT_ASSERT(0 != mesh);
-  *mesh = new Mesh(PETSC_COMM_WORLD, data.cellDim);
-  CPPUNIT_ASSERT(!mesh->isNull());
-  ALE::Obj<sieve_type> sieve = new sieve_type((*mesh)->comm());
+  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<ALE::Mesh::sieve_type> s = new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
-
+  ALE::Obj<ALE::Mesh::sieve_type> s = 
+    new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
+  
   ALE::SieveBuilder<ALE::Mesh>::buildTopology(s, data.cellDim, data.numCells,
-                                              const_cast<int*>(data.cells), data.numVertices,
+                                              const_cast<int*>(data.cells), 
+					      data.numVertices,
                                               interpolate, data.numCorners);
-  std::map<Mesh::point_type,Mesh::point_type> renumbering;
+  std::map<ALE::Mesh::point_type,ALE::Mesh::point_type> renumbering;
   ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
-  (*mesh)->setSieve(sieve);
-  (*mesh)->stratify();
-  ALE::SieveBuilder<Mesh>::buildCoordinates(*mesh, data.spaceDim, 
-					    data.vertices);
+  sieveMesh->setSieve(sieve);
+  sieveMesh->stratify();
+  ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, data.spaceDim, 
+						 data.vertices);
 
   // Material ids
-  const ALE::Obj<Mesh::label_sequence>& cells = (*mesh)->heightStratum(0);
+  const ALE::Obj<SieveMesh::label_sequence>& cells = 
+    sieveMesh->heightStratum(0);
   CPPUNIT_ASSERT(!cells.isNull());
-  const ALE::Obj<Mesh::label_type>& labelMaterials = 
-    (*mesh)->createLabel("material-id");
+  const ALE::Obj<SieveMesh::label_type>& labelMaterials = 
+    sieveMesh->createLabel("material-id");
   CPPUNIT_ASSERT(!labelMaterials.isNull());
   int i = 0;
-  for(Mesh::label_sequence::iterator e_iter=cells->begin(); 
+  for(SieveMesh::label_sequence::iterator e_iter=cells->begin(); 
       e_iter != cells->end();
       ++e_iter)
-    (*mesh)->setValue(labelMaterials, *e_iter, data.materialIds[i++]);
+    sieveMesh->setValue(labelMaterials, *e_iter, data.materialIds[i++]);
 
   // Groups
   for (int iGroup=0, index=0; iGroup < data.numGroups; ++iGroup) {
-    const ALE::Obj<int_section_type>& groupField = 
-      (*mesh)->getIntSection(data.groupNames[iGroup]);
+    const ALE::Obj<SieveMesh::int_section_type>& groupField = 
+      sieveMesh->getIntSection(data.groupNames[iGroup]);
     CPPUNIT_ASSERT(!groupField.isNull());
-    groupField->setChart((*mesh)->getSieve()->getChart());
+    groupField->setChart(sieveMesh->getSieve()->getChart());
 
     MeshIO::GroupPtType type;
     const int numPoints = data.groupSizes[iGroup];
@@ -89,14 +94,14 @@
         groupField->setFiberDimension(data.groups[index++], 1);
     } else if (0 == strcasecmp("vertex", data.groupTypes[iGroup])) {
       type = MeshIO::VERTEX;
-      const int numCells = (*mesh)->heightStratum(0)->size();
+      const int numCells = sieveMesh->heightStratum(0)->size();
       for(int i=0; i < numPoints; ++i)
         groupField->setFiberDimension(data.groups[index++]+numCells, 1);
     } else
       throw std::logic_error("Could not parse group type.");
-    (*mesh)->allocate(groupField);
+    sieveMesh->allocate(groupField);
   } // for
-  (*mesh)->getFactory()->clear();
+  sieveMesh->getFactory()->clear();
  
   return mesh;
 } // _createMesh
@@ -104,16 +109,22 @@
 // ----------------------------------------------------------------------
 // Check values in mesh against data.
 void
-pylith::meshio::TestMeshIO::_checkVals(const ALE::Obj<Mesh>& mesh,
+pylith::meshio::TestMeshIO::_checkVals(const topology::Mesh& mesh,
 				       const MeshData& data)
 { // _checkVals
+  typedef SieveMesh::int_section_type::chart_type chart_type;
+
   // Check mesh dimension
-  CPPUNIT_ASSERT_EQUAL(data.cellDim, mesh->getDimension());
+  CPPUNIT_ASSERT_EQUAL(data.cellDim, mesh.dimension());
 
+  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+  CPPUNIT_ASSERT(!sieveMesh.isNull());
+
   // Check vertices
-  const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
-  const ALE::Obj<Mesh::real_section_type>& coordsField =
-    mesh->getRealSection("coordinates");
+  const ALE::Obj<SieveMesh::label_sequence>& vertices = 
+    sieveMesh->depthStratum(0);
+  const ALE::Obj<SieveMesh::real_section_type>& coordsField =
+    sieveMesh->getRealSection("coordinates");
   const int numVertices = vertices->size();
   CPPUNIT_ASSERT(!vertices.isNull());
   CPPUNIT_ASSERT(!coordsField.isNull());
@@ -122,41 +133,41 @@
 		       coordsField->getFiberDimension(*vertices->begin()));
   int i = 0;
   const int spaceDim = data.spaceDim;
-  for(Mesh::label_sequence::iterator v_iter = 
+  for(SieveMesh::label_sequence::iterator v_iter = 
 	vertices->begin();
       v_iter != vertices->end();
       ++v_iter) {
-    const Mesh::real_section_type::value_type *vertexCoords = 
+    const SieveMesh::real_section_type::value_type *vertexCoords = 
       coordsField->restrictPoint(*v_iter);
     CPPUNIT_ASSERT(0 != vertexCoords);
     const double 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);
+				     tolerance);
       } else {
         CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vertexCoords[iDim]/data.vertices[i++],
-				   tolerance);
+				     tolerance);
       }
   } // for
 
   // check cells
-  const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
-  const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+  const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+  const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
 
   const int numCells = cells->size();
   CPPUNIT_ASSERT_EQUAL(data.numCells, numCells);
-  const int numCorners = mesh->getNumCellCorners();
+  const int numCorners = sieveMesh->getNumCellCorners();
   CPPUNIT_ASSERT_EQUAL(data.numCorners, numCorners);
 
-  ALE::ISieveVisitor::PointRetriever<Mesh::sieve_type> pV(sieve->getMaxConeSize());
+  ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> pV(sieve->getMaxConeSize());
   const int offset = numCells;
   i = 0;
-  for(Mesh::label_sequence::iterator e_iter = cells->begin();
+  for(SieveMesh::label_sequence::iterator e_iter = cells->begin();
       e_iter != cells->end();
       ++e_iter) {
     sieve->cone(*e_iter, pV);
-    const Mesh::point_type *cone = pV.getPoints();
+    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);
     }
@@ -164,23 +175,23 @@
   } // for
 
   // check materials
-  const ALE::Obj<Mesh::label_type>& labelMaterials = 
-    mesh->getLabel("material-id");
+  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(Mesh::label_sequence::iterator e_iter = cells->begin();
+  for(SieveMesh::label_sequence::iterator e_iter = cells->begin();
       e_iter != cells->end();
       ++e_iter)
-    materialIds[i++] = mesh->getValue(labelMaterials, *e_iter, idDefault);
+    materialIds[i++] = sieveMesh->getValue(labelMaterials, *e_iter, idDefault);
   
   for (int iCell=0; iCell < numCells; ++iCell)
     CPPUNIT_ASSERT_EQUAL(data.materialIds[iCell], materialIds[iCell]);
 
   // Check groups
   const ALE::Obj<std::set<std::string> >& groupNames = 
-    mesh->getIntSections();
+    sieveMesh->getIntSections();
   if (data.numGroups > 0) {
     CPPUNIT_ASSERT(!groupNames.isNull());
     CPPUNIT_ASSERT_EQUAL(data.numGroups, int(groupNames->size()));
@@ -190,22 +201,30 @@
   for (std::set<std::string>::const_iterator name=groupNames->begin();
        name != groupNames->end();
        ++name, ++iGroup) {
-    const ALE::Obj<int_section_type>& groupField = mesh->getIntSection(*name);
+    const ALE::Obj<SieveMesh::int_section_type>& groupField = 
+      sieveMesh->getIntSection(*name);
     CPPUNIT_ASSERT(!groupField.isNull());
-    const int_section_type::chart_type& chart = groupField->getChart();
-    Mesh::point_type firstPoint;
-    for(int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
-      if (groupField->getFiberDimension(*c_iter)) {firstPoint = *c_iter; break;}
-    }
+    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 = 
-      (mesh->height(firstPoint) == 0) ? "cell" : "vertex";
+      (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(int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
-      if (groupField->getFiberDimension(*c_iter)) points[i++] = *c_iter - offset;
-    }
+    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);

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIO.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIO.hh	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIO.hh	2009-01-24 00:22:09 UTC (rev 13936)
@@ -21,21 +21,25 @@
 #if !defined(pylith_meshio_testmeshio_hh)
 #define pylith_meshio_testmeshio_hh
 
+// Include directives ---------------------------------------------------
 #include <cppunit/extensions/HelperMacros.h>
 
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
-
+// Forward declarations -------------------------------------------------
 /// Namespace for pylith package
 namespace pylith {
   namespace meshio {
     class TestMeshIO;
-    class MeshIO;
+    class MeshIO; // USES MeshIO
 
-    class MeshData;
+    class MeshData; // test data
   } // meshio
+
+  namespace topology {
+    class Mesh; // USES Mesh
+  } // topology
 } // pylith
 
-/// C++ unit testing for TestMeshIO
+// MeshIO ---------------------------------------------------------------
 class pylith::meshio::TestMeshIO : public CppUnit::TestFixture
 { // class TestMeshIO
 
@@ -46,16 +50,16 @@
    *
    * @param data Mesh data
    *
-   * @returns PETSc mesh
+   * @returns PyLith mesh
    */
-  ALE::Obj<pylith::Mesh>* _createMesh(const MeshData& data);
+  topology::Mesh* _createMesh(const MeshData& data);
 
   /** Check values in mesh against data.
    *
-   * @param mesh PETSc mesh
+   * @param mesh PyLith mesh
    * @param data Mesh data
    */
-  void _checkVals(const ALE::Obj<Mesh>& mesh,
+  void _checkVals(const topology::Mesh& mesh,
 		  const MeshData& data);
 
   /** Test debug().
@@ -74,4 +78,5 @@
 
 #endif // pylith_meshio_testmeshio_hh
 
+
 // End of file 

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOAscii.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOAscii.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOAscii.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -16,7 +16,7 @@
 
 #include "pylith/meshio/MeshIOAscii.hh"
 
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include "pylith/topology/Mesh.hh" // USES Mesh
 
 #include "data/MeshData1D.hh"
 #include "data/MeshData1Din2D.hh"
@@ -155,7 +155,7 @@
 pylith::meshio::TestMeshIOAscii::_testWriteRead(const MeshData& data,
 						const char* filename)
 { // _testWriteRead
-  ALE::Obj<Mesh>* meshOut = _createMesh(data);
+  topology::Mesh* meshOut = _createMesh(data);
 
   // Write mesh
   MeshIOAscii iohandler;
@@ -164,7 +164,7 @@
   delete meshOut; meshOut = 0;
 
   // Read mesh
-  ALE::Obj<Mesh> meshIn;
+  topology::Mesh meshIn(PETSC_COMM_WORLD, data.cellDim);
   iohandler.read(&meshIn);
 
   // Make sure meshIn matches data
@@ -178,7 +178,7 @@
 					   const char* filename)
 { // _testWriteRead
   // Read mesh
-  ALE::Obj<Mesh> mesh;
+  topology::Mesh mesh(PETSC_COMM_WORLD, data.cellDim);
   MeshIOAscii iohandler;
   iohandler.filename(filename);
   iohandler.read(&mesh);

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOAscii.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOAscii.hh	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOAscii.hh	2009-01-24 00:22:09 UTC (rev 13936)
@@ -21,9 +21,10 @@
 #if !defined(pylith_meshio_testmeshioascii_hh)
 #define pylith_meshio_testmeshioascii_hh
 
+// Include directives ---------------------------------------------------
 #include "TestMeshIO.hh"
 
-/// Namespace for pylith package
+// Forward declarations -------------------------------------------------
 namespace pylith {
   namespace meshio {
     class TestMeshIOAscii;
@@ -31,7 +32,7 @@
   } // meshio
 } // pylith
 
-/// C++ unit testing for Quadrature1D
+// TestMeshIOAscii ------------------------------------------------------
 class pylith::meshio::TestMeshIOAscii : public TestMeshIO
 { // class TestMeshIOAscii
 

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOCubit.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOCubit.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOCubit.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -16,7 +16,7 @@
 
 #include "pylith/meshio/MeshIOCubit.hh"
 
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/utils/array.hh" // USES int_array
 
 #include "data/MeshDataCubitTri.hh"
@@ -117,7 +117,7 @@
   iohandler.filename(filename);
 
   // Read mesh
-  ALE::Obj<Mesh> mesh;
+  topology::Mesh mesh(PETSC_COMM_WORLD, data.cellDim);
   iohandler.read(&mesh);
 
   // Make sure meshIn matches data

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOCubit.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOCubit.hh	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOCubit.hh	2009-01-24 00:22:09 UTC (rev 13936)
@@ -21,9 +21,10 @@
 #if !defined(pylith_meshio_testmeshiocubit_hh)
 #define pylith_meshio_testmeshiocubit_hh
 
+// Include directives ---------------------------------------------------
 #include "TestMeshIO.hh"
 
-/// Namespace for pylith package
+// Forward declarations -------------------------------------------------
 namespace pylith {
   namespace meshio {
     class TestMeshIOCubit;
@@ -31,7 +32,7 @@
   } // meshio
 } // pylith
 
-/// C++ unit testing for Quadrature1D
+// TestMeshIOCubit ------------------------------------------------------
 class pylith::meshio::TestMeshIOCubit : public TestMeshIO
 { // class TestMeshIOCubit
 

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOLagrit.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOLagrit.cc	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOLagrit.cc	2009-01-24 00:22:09 UTC (rev 13936)
@@ -16,7 +16,7 @@
 
 #include "pylith/meshio/MeshIOLagrit.hh"
 
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/utils/array.hh" // USES int_array
 
 #include "data/MeshDataLagritTet.hh"
@@ -146,7 +146,7 @@
 #endif
 
   // Read mesh
-  ALE::Obj<Mesh> mesh;
+  topology::Mesh mesh(PETSC_COMM_WORLD, data.cellDim);
   iohandler.read(&mesh);
 
   // Make sure meshIn matches data

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOLagrit.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOLagrit.hh	2009-01-24 00:20:23 UTC (rev 13935)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/meshio/TestMeshIOLagrit.hh	2009-01-24 00:22:09 UTC (rev 13936)
@@ -21,9 +21,10 @@
 #if !defined(pylith_meshio_testmeshiolagrit_hh)
 #define pylith_meshio_testmeshiolagrit_hh
 
+// Include directives ---------------------------------------------------
 #include "TestMeshIO.hh"
 
-/// Namespace for pylith package
+// Forward declarations -------------------------------------------------
 namespace pylith {
   namespace meshio {
     class TestMeshIOLagrit;
@@ -31,7 +32,7 @@
   } // meshio
 } // pylith
 
-/// C++ unit testing for Quadrature1D
+// TestMeshIOLagrit -----------------------------------------------------
 class pylith::meshio::TestMeshIOLagrit : public TestMeshIO
 { // class TestMeshIOLagrit
 



More information about the CIG-COMMITS mailing list