[cig-commits] r21720 - short/3D/PyLith/trunk/unittests/libtests/faults

brad at geodynamics.org brad at geodynamics.org
Thu Apr 4 11:29:47 PDT 2013


Author: brad
Date: 2013-04-04 11:29:46 -0700 (Thu, 04 Apr 2013)
New Revision: 21720

Added:
   short/3D/PyLith/trunk/unittests/libtests/faults/TestSlipFn.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestSlipFn.hh
Modified:
   short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
   short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.hh
   short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.hh
   short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.hh
   short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.hh
   short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.hh
Log:
Code cleanup of slip time function tests. Updated to use visitors.

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am	2013-04-04 18:29:46 UTC (rev 21720)
@@ -27,14 +27,15 @@
 
 # Primary source files
 testfaults_SOURCES = \
-	TestFault.cc \
-	TestFaultCohesive.cc \
+	TestSlipFn.cc \
 	TestStepSlipFn.cc \
 	TestConstRateSlipFn.cc \
 	TestBruneSlipFn.cc \
 	TestLiuCosSlipFn.cc \
 	TestTimeHistorySlipFn.cc \
 	TestEqKinSrc.cc \
+	TestFault.cc \
+	TestFaultCohesive.cc \
 	TestFaultCohesiveKin.cc \
 	TestFaultCohesiveKinLine2.cc \
 	TestFaultCohesiveKinTri3.cc \
@@ -64,6 +65,7 @@
 
 
 noinst_HEADERS = \
+	TestSlipFn.hh \
 	TestBruneSlipFn.hh \
 	TestLiuCosSlipFn.hh \
 	TestTimeHistorySlipFn.hh \

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc	2013-04-04 18:29:46 UTC (rev 21720)
@@ -22,7 +22,6 @@
 
 #include "pylith/faults/BruneSlipFn.hh" // USES BruneSlipFn
 
-#include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Field.hh" // USES Field
@@ -40,11 +39,6 @@
 CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestBruneSlipFn );
 
 // ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
 namespace pylith {
   namespace faults {
     namespace _TestBruneSlipFn {
@@ -70,7 +64,11 @@
 void
 pylith::faults::TestBruneSlipFn::testConstructor(void)
 { // testConstructor
+  PYLITH_METHOD_BEGIN;
+
   BruneSlipFn slipfn;
+
+  PYLITH_METHOD_END;
 } // testConstructor
 
 // ----------------------------------------------------------------------
@@ -78,6 +76,8 @@
 void
 pylith::faults::TestBruneSlipFn::testDbFinalSlip(void)
 { // testDbFinalSlip
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABC";
   BruneSlipFn slipfn;
   
@@ -88,6 +88,8 @@
   CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbFinalSlip->label()));
   CPPUNIT_ASSERT(!slipfn._dbSlipTime);
   CPPUNIT_ASSERT(!slipfn._dbRiseTime);
+
+  PYLITH_METHOD_END;
 } // testDbFinalSlip
 
 // ----------------------------------------------------------------------
@@ -95,6 +97,8 @@
 void
 pylith::faults::TestBruneSlipFn::testDbSlipTime(void)
 { // testDbSlipTime
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABCD";
   BruneSlipFn slipfn;
   
@@ -105,6 +109,8 @@
   CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbSlipTime->label()));
   CPPUNIT_ASSERT(!slipfn._dbFinalSlip);
   CPPUNIT_ASSERT(!slipfn._dbRiseTime);
+
+  PYLITH_METHOD_END;
 } // testDbSlipTime
 
 // ----------------------------------------------------------------------
@@ -112,6 +118,8 @@
 void
 pylith::faults::TestBruneSlipFn::testDbRiseTime(void)
 { // testDbRiseTime
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABCDE";
   BruneSlipFn slipfn;
   
@@ -119,10 +127,11 @@
   slipfn.dbRiseTime(&db);
 
   CPPUNIT_ASSERT(slipfn._dbRiseTime);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbRiseTime->label()));
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbRiseTime->label()));
   CPPUNIT_ASSERT(!slipfn._dbFinalSlip);
   CPPUNIT_ASSERT(!slipfn._dbSlipTime);
+
+  PYLITH_METHOD_END;
 } // testDbRiseTime
 
 // ----------------------------------------------------------------------
@@ -130,6 +139,8 @@
 void
 pylith::faults::TestBruneSlipFn::testInitialize1D(void)
 { // testInitialize1D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/line2.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -154,6 +165,8 @@
 				       riseTimeE,
 				       numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize1D
 
 // ----------------------------------------------------------------------
@@ -161,6 +174,8 @@
 void
 pylith::faults::TestBruneSlipFn::testInitialize2D(void)
 { // testInitialize2D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tri3.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -186,6 +201,8 @@
 				       riseTimeE,
 				       numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize2D
 
 // ----------------------------------------------------------------------
@@ -193,6 +210,8 @@
 void
 pylith::faults::TestBruneSlipFn::testInitialize3D(void)
 { // testInitialize3D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tet4.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -219,6 +238,8 @@
 				       riseTimeE,
 				       numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize3D
 
 // ----------------------------------------------------------------------
@@ -226,6 +247,8 @@
 void
 pylith::faults::TestBruneSlipFn::testSlip(void)
 { // testSlip
+  PYLITH_METHOD_BEGIN;
+
   const PylithScalar finalSlipE[] = { 2.3, 0.1, 
 				0.0, 0.0};
   const PylithScalar slipTimeE[] = { 1.2, 1.3 };
@@ -258,8 +281,9 @@
   const PylithScalar tolerance = 1.0e-06;
   for(PetscInt v = vStart, iPoint=0; v < vEnd; ++v, ++iPoint) {
     PylithScalar slipMag = 0.0;
-    for (int iDim=0; iDim < spaceDim; ++iDim)
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
       slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
+    } // for
     slipMag = sqrt(slipMag);
     const PylithScalar peakRate = slipMag / riseTimeE[iPoint] * 1.745;
     const PylithScalar tau = (slipMag > 0.0) ? slipMag / (exp(1.0) * peakRate) : 1.0;
@@ -274,6 +298,8 @@
       CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, slipArray[off+d], tolerance);
     } // for
   } // for
+
+  PYLITH_METHOD_END;
 } // testSlip
 
 // ----------------------------------------------------------------------
@@ -281,6 +307,8 @@
 void
 pylith::faults::TestBruneSlipFn::testSlipTH(void)
 { // testSlipTH
+  PYLITH_METHOD_BEGIN;
+
   const PylithScalar t = 0.734;
   const PylithScalar finalSlip = 4.64;
   const PylithScalar riseTime = 3.23;
@@ -299,6 +327,8 @@
 
   slip = BruneSlipFn::_slipFn(1.0e+10, finalSlip, riseTime);
   CPPUNIT_ASSERT_DOUBLES_EQUAL(finalSlip, slip, tolerance);
+
+  PYLITH_METHOD_END;
 } // testSlipTH
 
 // ----------------------------------------------------------------------
@@ -309,6 +339,8 @@
 					     BruneSlipFn* slipfn,
 					     const PylithScalar originTime)
 { // _initialize
+  PYLITH_METHOD_BEGIN;
+
   CPPUNIT_ASSERT(mesh);
   CPPUNIT_ASSERT(faultMesh);
   CPPUNIT_ASSERT(slipfn);
@@ -335,33 +367,8 @@
   mesh->coordsys(&cs);
 
   // Create fault mesh
-  PetscDM dmMesh = mesh->dmMesh(), faultBoundaryDM;
-  PetscInt firstFaultVertex = 0;
-  PetscInt firstLagrangeVertex, firstFaultCell;
-  DMLabel groupField;
-  const bool useLagrangeConstraints = true;
+  TestSlipFn::_createFaultMesh(faultMesh, mesh, faultLabel, faultId);
 
-  err = DMPlexGetStratumSize(dmMesh, faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
-  firstFaultCell = firstLagrangeVertex;
-  if (useLagrangeConstraints) {
-    firstFaultCell += firstLagrangeVertex;
-  } // if
-  err = DMPlexGetLabel(dmMesh, faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  CohesiveTopology::createFault(faultMesh, faultBoundary, faultBoundaryDM,
-                                *mesh, groupField);
-  CohesiveTopology::create(mesh, *faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           faultId,
-                           firstFaultVertex, firstLagrangeVertex, firstFaultCell,
-                           useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are not
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(mesh, faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
   spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -386,6 +393,8 @@
   slipfn->dbRiseTime(&dbRiseTime);
   
   slipfn->initialize(*faultMesh, normalizer, originTime);
+
+  PYLITH_METHOD_END;
 } // _initialize
 
 // ----------------------------------------------------------------------
@@ -393,6 +402,8 @@
 void
 pylith::faults::TestBruneSlipFn::_testInitialize(const _TestBruneSlipFn::DataStruct& data)
 { // _testInitialize
+  PYLITH_METHOD_BEGIN;
+
   PetscErrorCode err;
   // Setup mesh
   topology::Mesh mesh;
@@ -411,33 +422,8 @@
 
   // Create fault mesh
   topology::SubMesh faultMesh;
-  DM       dmMesh = mesh.dmMesh(), faultBoundaryDM;
-  PetscInt firstFaultVertex = 0;
-  PetscInt firstLagrangeVertex, firstFaultCell;
-  DMLabel  groupField;
-  const bool useLagrangeConstraints = true;
+  TestSlipFn::_createFaultMesh(&faultMesh, &mesh, data.faultLabel, data.faultId);
 
-  err = DMPlexGetStratumSize(dmMesh, data.faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
-  firstFaultCell = firstLagrangeVertex;
-  if (useLagrangeConstraints) {
-    firstFaultCell += firstLagrangeVertex;
-  }
-  err = DMPlexGetLabel(dmMesh, data.faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  CohesiveTopology::createFault(&faultMesh, faultBoundary, faultBoundaryDM,
-                                mesh, groupField);
-  CohesiveTopology::create(&mesh, faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           data.faultId,
-                           firstFaultVertex, firstLagrangeVertex, firstFaultCell,
-                           useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are not
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(&mesh, &faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
   spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -466,109 +452,40 @@
   slipfn.initialize(faultMesh, normalizer, originTime);
 
   CPPUNIT_ASSERT(slipfn._parameters);
-  PetscSection finalSlipSection = slipfn._parameters->get("final slip").petscSection();
-  Vec          finalSlipVec     = slipfn._parameters->get("final slip").localVector();
-  PetscScalar *finalSlipArray;
-  PetscSection slipTimeSection = slipfn._parameters->get("slip time").petscSection();
-  Vec          slipTimeVec     = slipfn._parameters->get("slip time").localVector();
-  PetscScalar *slipTimeArray;
-  CPPUNIT_ASSERT(slipTimeSection);CPPUNIT_ASSERT(slipTimeVec);
-  PetscSection riseTimeSection = slipfn._parameters->get("rise time").petscSection();
-  Vec          riseTimeVec     = slipfn._parameters->get("rise time").localVector();
-  PetscScalar *riseTimeArray;
-  CPPUNIT_ASSERT(riseTimeSection);CPPUNIT_ASSERT(riseTimeVec);
+  topology::VecVisitorMesh finalSlipVisitor(slipfn._parameters->get("final slip"));
+  const PetscScalar* finalSlipArray = finalSlipVisitor.localArray();CPPUNIT_ASSERT(finalSlipArray);
 
-  const PylithScalar tolerance = 1.0e-06;
-  PetscInt           vStart, vEnd, iPoint = 0;
-  err = DMPlexGetDepthStratum(faultMesh.dmMesh(), 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(slipTimeVec,  &slipTimeArray);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(riseTimeVec,  &riseTimeArray);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
-    PetscInt fsdof, fsoff, stdof, stoff, rtdof, rtoff;
+  topology::VecVisitorMesh slipTimeVisitor(slipfn._parameters->get("slip time"));
+  const PetscScalar* slipTimeArray = slipTimeVisitor.localArray();CPPUNIT_ASSERT(slipTimeArray);
 
-    err = PetscSectionGetDof(finalSlipSection, v, &fsdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(finalSlipSection, v, &fsoff);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetDof(riseTimeSection, v, &rtdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(riseTimeSection, v, &rtoff);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(spaceDim, fsdof);
-    for(PetscInt d = 0; d < fsdof; ++d)
-      CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+d], finalSlipArray[fsoff+d], tolerance);
+  topology::VecVisitorMesh riseTimeVisitor(slipfn._parameters->get("rise time"));
+  const PetscScalar* riseTimeArray = riseTimeVisitor.localArray();CPPUNIT_ASSERT(riseTimeArray);
 
-    CPPUNIT_ASSERT_EQUAL(1, stdof);
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
+  PetscDM faultDMMesh = faultMesh.dmMesh();CPPUNIT_ASSERT(faultDMMesh);
+  topology::Stratum verticesStratum(faultDMMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
 
-    CPPUNIT_ASSERT_EQUAL(1, rtdof);
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(data.riseTimeE[iPoint], riseTimeArray[rtoff], tolerance);
-  } // for
-  err = VecRestoreArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(slipTimeVec,  &slipTimeArray);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(riseTimeVec,  &riseTimeArray);CHECK_PETSC_ERROR(err);
-} // _testInitialize
+  const PylithScalar tolerance = 1.0e-06;
+  for(PetscInt v = vStart, iPoint = 0; v < vEnd; ++v, ++iPoint) {
+    const PetscInt fsoff = finalSlipVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, finalSlipVisitor.sectionDof(v));
 
-// ----------------------------------------------------------------------
-// Setup fault coordinates
-void
-pylith::faults::TestBruneSlipFn::_setupFaultCoordinates(topology::Mesh *mesh, topology::SubMesh *faultMesh)
-{ // _setupFaultCoordinates
-  const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
-  if (!faultSieveMesh.isNull()) {
-    const ALE::Obj<RealSection>& oldCoordSection = mesh->sieveMesh()->getRealSection("coordinates");
-    faultSieveMesh->setRealSection("coordinates", oldCoordSection);
-  }
+    const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(1, slipTimeVisitor.sectionDof(v));
 
-  DM              dmMesh      = mesh->dmMesh();
-  DM              faultDMMesh = faultMesh->dmMesh();
-  const PetscInt  spaceDim    = mesh->dimension();
-  IS              subpointIS;
-  const PetscInt *points;
-  PetscSection    coordSection, fcoordSection;
-  PetscInt        vStart, vEnd, ffStart, ffEnd;
-  PetscErrorCode  err;
+    const PetscInt rtoff = riseTimeVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(1, riseTimeVisitor.sectionDof(v));
 
-  CPPUNIT_ASSERT(dmMesh);
-  CPPUNIT_ASSERT(faultDMMesh);
-  err = DMPlexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetHeightStratum(faultDMMesh, 1, &ffStart, &ffEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexCreateSubpointIS(faultDMMesh, &subpointIS);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(faultDMMesh, &fcoordSection);CHECK_PETSC_ERROR(err);
-  err = PetscSectionSetChart(fcoordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    err = PetscSectionSetDof(fcoordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
-  }
-  err = PetscSectionSetUp(fcoordSection);CHECK_PETSC_ERROR(err);
-  Vec          coordVec, fcoordVec;
-  PetscScalar *coords,  *fcoords;
-  PetscInt     coordSize;
-
-  err = PetscSectionGetStorageSize(fcoordSection, &coordSize);CHECK_PETSC_ERROR(err);
-  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
-  err = VecCreate(mesh->comm(), &fcoordVec);CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(fcoordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(fcoordVec);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    PetscInt off, foff;
-
-    // Notice that subpointMap[] does not account for cohesive cells
-    err = PetscSectionGetOffset(coordSection, points[v]+(ffEnd-ffStart), &off);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(fcoordSection, v, &foff);CHECK_PETSC_ERROR(err);
     for(PetscInt d = 0; d < spaceDim; ++d) {
-      fcoords[foff+d] = coords[off+d];
-    }
-  }
-  err = ISRestoreIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  err = DMSetCoordinatesLocal(faultDMMesh, fcoordVec);CHECK_PETSC_ERROR(err);
-} // _setupFaultCoordinates
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+d], finalSlipArray[fsoff+d], tolerance);
+    } // for
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(data.riseTimeE[iPoint], riseTimeArray[rtoff], tolerance);
+  } // for
 
+  PYLITH_METHOD_END;
+} // _testInitialize
 
 
 // End of file 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.hh	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.hh	2013-04-04 18:29:46 UTC (rev 21720)
@@ -27,11 +27,11 @@
 #if !defined(pylith_faults_testbruneslipfn_hh)
 #define pylith_faults_testbruneslipfn_hh
 
+#include "TestSlipFn.hh" // ISA TestSlipFn
+
 #include "pylith/faults/faultsfwd.hh" // USES BruneSlipFn
 #include "pylith/topology/topologyfwd.hh" // USES Mesh, SubMesh
 
-#include <cppunit/extensions/HelperMacros.h>
-
 /// Namespace for pylith package
 namespace pylith {
   namespace faults {
@@ -44,7 +44,7 @@
 } // pylith
 
 /// C++ unit testing for BruneSlipFn
-class pylith::faults::TestBruneSlipFn : public CppUnit::TestFixture
+class pylith::faults::TestBruneSlipFn : public TestSlipFn
 { // class TestBruneSlipFn
 
   // CPPUNIT TEST SUITE /////////////////////////////////////////////////

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc	2013-04-04 18:29:46 UTC (rev 21720)
@@ -22,11 +22,12 @@
 
 #include "pylith/faults/ConstRateSlipFn.hh" // USES ConstRateSlipFn
 
-#include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 
 #include "spatialdata/geocoords/CSCart.hh" // USES CSCart
@@ -57,16 +58,15 @@
 } // pylith
 
 // ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
 // Test constructor.
 void
 pylith::faults::TestConstRateSlipFn::testConstructor(void)
 { // testConstructor
+  PYLITH_METHOD_BEGIN;
+
   ConstRateSlipFn slipfn;
+
+  PYLITH_METHOD_END;
 } // testConstructor
 
 // ----------------------------------------------------------------------
@@ -74,16 +74,19 @@
 void
 pylith::faults::TestConstRateSlipFn::testDbSlipRate(void)
 { // testDbSlipRate
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABC";
   ConstRateSlipFn slipfn;
   
   spatialdata::spatialdb::SimpleDB db(label);
   slipfn.dbSlipRate(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbSlipRate);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbSlipRate->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbSlipTime);
+  CPPUNIT_ASSERT(slipfn._dbSlipRate);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbSlipRate->label()));
+  CPPUNIT_ASSERT(!slipfn._dbSlipTime);
+
+  PYLITH_METHOD_END;
 } // testDbSlipRate
 
 // ----------------------------------------------------------------------
@@ -91,16 +94,19 @@
 void
 pylith::faults::TestConstRateSlipFn::testDbSlipTime(void)
 { // testDbSlipTime
-  const char* label = "database ABCD";
+  PYLITH_METHOD_BEGIN;
+
+ const char* label = "database ABCD";
   ConstRateSlipFn slipfn;
   
   spatialdata::spatialdb::SimpleDB db(label);
   slipfn.dbSlipTime(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbSlipTime);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbSlipTime->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbSlipRate);
+  CPPUNIT_ASSERT(slipfn._dbSlipTime);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbSlipTime->label()));
+  CPPUNIT_ASSERT(!slipfn._dbSlipRate);
+
+  PYLITH_METHOD_END;
 } // testDbSlipTime
 
 // ----------------------------------------------------------------------
@@ -108,6 +114,8 @@
 void
 pylith::faults::TestConstRateSlipFn::testInitialize1D(void)
 { // testInitialize1D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/line2.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -128,6 +136,8 @@
 					   slipTimeE,
 					   numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize1D
 
 // ----------------------------------------------------------------------
@@ -135,6 +145,8 @@
 void
 pylith::faults::TestConstRateSlipFn::testInitialize2D(void)
 { // testInitialize2D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tri3.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -156,6 +168,8 @@
 					   slipTimeE,
 					   numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize2D
 
 // ----------------------------------------------------------------------
@@ -163,6 +177,8 @@
 void
 pylith::faults::TestConstRateSlipFn::testInitialize3D(void)
 { // testInitialize3D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tet4.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -185,6 +201,8 @@
 					   slipTimeE,
 					   numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize3D
 
 // ----------------------------------------------------------------------
@@ -192,6 +210,8 @@
 void
 pylith::faults::TestConstRateSlipFn::testSlip(void)
 { // testSlip
+  PYLITH_METHOD_BEGIN;
+
   const PylithScalar slipRateE[] = { 0.1, 0.2, 
 			       0.3, 0.4};
   const PylithScalar slipTimeE[] = { 1.2, 1.3 };
@@ -202,16 +222,9 @@
   ConstRateSlipFn slipfn;
   _initialize(&mesh, &faultMesh, &slipfn, originTime);
   
-  const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
-  CPPUNIT_ASSERT(0 != cs);
-
+  const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();CPPUNIT_ASSERT(cs);
   const int spaceDim = cs->spaceDim();
-  DM dmMesh = faultMesh.dmMesh();
-  PetscErrorCode err;
-  CPPUNIT_ASSERT(dmMesh);
 
-  PetscInt       vStart, vEnd;
-  err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
   topology::Field<topology::SubMesh> slip(faultMesh);
   slip.newSection(topology::Field<topology::Mesh>::VERTICES_FIELD, spaceDim);
   slip.allocate();
@@ -219,26 +232,28 @@
   const PylithScalar t = 1.234;
   slipfn.slip(&slip, originTime+t);
 
+  PetscDM dmMesh = faultMesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
+
+  topology::VecVisitorMesh slipVisitor(slip);
+  const PetscScalar* slipArray = slipVisitor.localArray();CPPUNIT_ASSERT(slipArray);
+
   const PylithScalar tolerance = 1.0e-06;
-  int iPoint = 0;
-  PetscSection slipSection = slip.petscSection();
-  Vec          slipVec     = slip.localVector();
-  PetscScalar *slipArray;
-  CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
-  err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
+  for(PetscInt v = vStart, iPoint = 0; v < vEnd; ++v, ++iPoint) {
     const PylithScalar t0 = slipTimeE[iPoint];
-    PetscInt dof, off;
 
-    err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
-    for(PetscInt d = 0; d < dof; ++d) {
+    const PetscInt off = slipVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, slipVisitor.sectionDof(v));
+
+    for(PetscInt d = 0; d < spaceDim; ++d) {
       const PylithScalar slipE = (t - t0) > 0.0 ? slipRateE[iPoint*spaceDim+d] * (t - t0) : 0.0;
       CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, slipArray[off+d], tolerance);
     }
   } // for
-  err = VecRestoreArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+
+  PYLITH_METHOD_END;
 } // testSlip
 
 // ----------------------------------------------------------------------
@@ -249,9 +264,11 @@
 					    	 ConstRateSlipFn* slipfn,
 						 const PylithScalar originTime)
 { // _initialize
-  CPPUNIT_ASSERT(0 != mesh);
-  CPPUNIT_ASSERT(0 != faultMesh);
-  CPPUNIT_ASSERT(0 != slipfn);
+  PYLITH_METHOD_BEGIN;
+
+  CPPUNIT_ASSERT(mesh);
+  CPPUNIT_ASSERT(faultMesh);
+  CPPUNIT_ASSERT(slipfn);
   PetscErrorCode err;
 
   const char* meshFilename = "data/tri3.mesh";
@@ -274,33 +291,8 @@
   mesh->coordsys(&cs);
 
   // Create fault mesh
-  DM       dmMesh = mesh->dmMesh(), faultBoundaryDM;
-  PetscInt firstFaultVertex = 0;
-  PetscInt firstLagrangeVertex, firstFaultCell;
-  DMLabel  groupField;
-  const bool useLagrangeConstraints = true;
+  TestSlipFn::_createFaultMesh(faultMesh, mesh, faultLabel, faultId);
 
-  err = DMPlexGetStratumSize(dmMesh, faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
-  firstFaultCell = firstLagrangeVertex;
-  if (useLagrangeConstraints) {
-    firstFaultCell += firstLagrangeVertex;
-  }
-  err = DMPlexGetLabel(dmMesh, faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  CohesiveTopology::createFault(faultMesh, faultBoundary, faultBoundaryDM,
-                                *mesh, groupField);
-  CohesiveTopology::create(mesh, *faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           faultId,
-                           firstFaultVertex, firstLagrangeVertex, firstFaultCell,
-                           useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are not
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(mesh, faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbSlipRate("slip rate");
   spatialdata::spatialdb::SimpleIOAscii ioSlipRate;
@@ -319,6 +311,8 @@
   slipfn->dbSlipTime(&dbSlipTime);
   
   slipfn->initialize(*faultMesh, normalizer, originTime);
+
+  PYLITH_METHOD_END;
 } // _initialize
 
 // ----------------------------------------------------------------------
@@ -326,6 +320,8 @@
 void
 pylith::faults::TestConstRateSlipFn::_testInitialize(const _TestConstRateSlipFn::DataStruct& data)
 { // _testInitialize
+  PYLITH_METHOD_BEGIN;
+
   PetscErrorCode err;
   // Setup mesh
   topology::Mesh mesh;
@@ -344,33 +340,8 @@
 
   // Create fault mesh
   topology::SubMesh faultMesh;
-  DM       dmMesh = mesh.dmMesh(), faultBoundaryDM;
-  PetscInt firstFaultVertex = 0;
-  PetscInt firstLagrangeVertex, firstFaultCell;
-  DMLabel  groupField;
-  const bool useLagrangeConstraints = true;
+  TestSlipFn::_createFaultMesh(&faultMesh, &mesh, data.faultLabel, data.faultId);
 
-  err = DMPlexGetStratumSize(dmMesh, data.faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
-  firstFaultCell = firstLagrangeVertex;
-  if (useLagrangeConstraints) {
-    firstFaultCell += firstLagrangeVertex;
-  }
-  err = DMPlexGetLabel(dmMesh, data.faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  CohesiveTopology::createFault(&faultMesh, faultBoundary, faultBoundaryDM,
-                                mesh, groupField);
-  CohesiveTopology::create(&mesh, faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           data.faultId,
-                           firstFaultVertex, firstLagrangeVertex, firstFaultCell,
-                           useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are not
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(&mesh, &faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbSlipRate("slip rate");
   spatialdata::spatialdb::SimpleIOAscii ioSlipRate;
@@ -392,100 +363,34 @@
   
   slipfn.initialize(faultMesh, normalizer, originTime);
 
-  CPPUNIT_ASSERT(0 != slipfn._parameters);
-  PetscSection slipRateSection = slipfn._parameters->get("slip rate").petscSection();
-  Vec          slipRateVec     = slipfn._parameters->get("slip rate").localVector();
-  PetscScalar *slipRateArray;
-  CPPUNIT_ASSERT(slipRateSection);CPPUNIT_ASSERT(slipRateVec);
-  PetscSection slipTimeSection = slipfn._parameters->get("slip time").petscSection();
-  Vec          slipTimeVec     = slipfn._parameters->get("slip time").localVector();
-  PetscScalar *slipTimeArray;
-  CPPUNIT_ASSERT(slipTimeSection);CPPUNIT_ASSERT(slipTimeVec);
+  CPPUNIT_ASSERT(slipfn._parameters);
+  topology::VecVisitorMesh slipRateVisitor(slipfn._parameters->get("slip rate"));
+  const PetscScalar* slipRateArray = slipRateVisitor.localArray();CPPUNIT_ASSERT(slipRateArray);
 
+  topology::VecVisitorMesh slipTimeVisitor(slipfn._parameters->get("slip time"));
+  const PetscScalar* slipTimeArray = slipTimeVisitor.localArray();CPPUNIT_ASSERT(slipTimeArray);
+
+  PetscDM faultDMMesh = faultMesh.dmMesh();CPPUNIT_ASSERT(faultDMMesh);
+  topology::Stratum verticesStratum(faultDMMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
+
   const PylithScalar tolerance = 1.0e-06;
-  PetscInt           vStart, vEnd, iPoint = 0;
-  err = DMPlexGetDepthStratum(faultMesh.dmMesh(), 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(slipRateVec, &slipRateArray);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
-    PetscInt srdof, sroff, stdof, stoff;
+  for(PetscInt v = vStart, iPoint = 0; v < vEnd; ++v, ++iPoint) {
+    const PetscInt sroff = slipRateVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, slipRateVisitor.sectionDof(v));
 
-    err = PetscSectionGetDof(slipRateSection, v, &srdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipRateSection, v, &sroff);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(spaceDim, srdof);
+    const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(1, slipTimeVisitor.sectionDof(v));
+
     for(PetscInt d = 0; d < spaceDim; ++d)
       CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipRateE[iPoint*spaceDim+d], slipRateArray[sroff+d], tolerance);
 
-    CPPUNIT_ASSERT_EQUAL(1, stdof);
     CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
   } // for
-  err = VecRestoreArray(slipRateVec, &slipRateArray);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
+
+  PYLITH_METHOD_END;
 } // _testInitialize
 
-// ----------------------------------------------------------------------
-// Setup fault coordinates
-void
-pylith::faults::TestConstRateSlipFn::_setupFaultCoordinates(topology::Mesh *mesh, topology::SubMesh *faultMesh)
-{ // _setupFaultCoordinates
-  const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
-  if (!faultSieveMesh.isNull()) {
-    const ALE::Obj<RealSection>& oldCoordSection = mesh->sieveMesh()->getRealSection("coordinates");
-    faultSieveMesh->setRealSection("coordinates", oldCoordSection);
-  }
 
-  DM              dmMesh      = mesh->dmMesh();
-  DM              faultDMMesh = faultMesh->dmMesh();
-  const PetscInt  spaceDim    = mesh->dimension();
-  IS              subpointIS;
-  const PetscInt *points;
-  PetscSection    coordSection, fcoordSection;
-  PetscInt        vStart, vEnd, ffStart, ffEnd;
-  PetscErrorCode  err;
-
-  CPPUNIT_ASSERT(dmMesh);
-  CPPUNIT_ASSERT(faultDMMesh);
-  err = DMPlexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetHeightStratum(faultDMMesh, 1, &ffStart, &ffEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexCreateSubpointIS(faultDMMesh, &subpointIS);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(faultDMMesh, &fcoordSection);CHECK_PETSC_ERROR(err);
-  err = PetscSectionSetChart(fcoordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    err = PetscSectionSetDof(fcoordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
-  }
-  err = PetscSectionSetUp(fcoordSection);CHECK_PETSC_ERROR(err);
-  Vec          coordVec, fcoordVec;
-  PetscScalar *coords,  *fcoords;
-  PetscInt     coordSize;
-
-  err = PetscSectionGetStorageSize(fcoordSection, &coordSize);CHECK_PETSC_ERROR(err);
-  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
-  err = VecCreate(mesh->comm(), &fcoordVec);CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(fcoordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(fcoordVec);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    PetscInt off, foff;
-
-    // Notice that subpointMap[] does not account for cohesive cells
-    err = PetscSectionGetOffset(coordSection, points[v]+(ffEnd-ffStart), &off);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(fcoordSection, v, &foff);CHECK_PETSC_ERROR(err);
-    for(PetscInt d = 0; d < spaceDim; ++d) {
-      fcoords[foff+d] = coords[off+d];
-    }
-  }
-  err = ISRestoreIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  err = DMSetCoordinatesLocal(faultDMMesh, fcoordVec);CHECK_PETSC_ERROR(err);
-} // _setupFaultCoordinates
-
-
-
 // End of file 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.hh	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.hh	2013-04-04 18:29:46 UTC (rev 21720)
@@ -27,11 +27,11 @@
 #if !defined(pylith_faults_testconstrateslipfn_hh)
 #define pylith_faults_testconstrateslipfn_hh
 
+#include "TestSlipFn.hh" // ISA TestSlipFn
+
 #include "pylith/faults/faultsfwd.hh" // USES ConstRateSlipFn
 #include "pylith/topology/topologyfwd.hh" // USES Mesh, SubMesh
 
-#include <cppunit/extensions/HelperMacros.h>
-
 /// Namespace for pylith package
 namespace pylith {
   namespace faults {
@@ -44,7 +44,7 @@
 } // pylith
 
 /// C++ unit testing for ConstRateSlipFn
-class pylith::faults::TestConstRateSlipFn : public CppUnit::TestFixture
+class pylith::faults::TestConstRateSlipFn : public TestSlipFn
 { // class TestConstRateSlipFn
 
   // CPPUNIT TEST SUITE /////////////////////////////////////////////////

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc	2013-04-04 18:29:46 UTC (rev 21720)
@@ -22,11 +22,12 @@
 
 #include "pylith/faults/LiuCosSlipFn.hh" // USES LiuCosSlipFn
 
-#include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 
 #include "spatialdata/geocoords/CSCart.hh" // USES CSCart
@@ -38,11 +39,6 @@
 CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestLiuCosSlipFn );
 
 // ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
 namespace pylith {
   namespace faults {
     namespace _TestLiuCosSlipFn {
@@ -68,7 +64,11 @@
 void
 pylith::faults::TestLiuCosSlipFn::testConstructor(void)
 { // testConstructor
+  PYLITH_METHOD_BEGIN;
+
   LiuCosSlipFn slipfn;
+
+  PYLITH_METHOD_END;
 } // testConstructor
 
 // ----------------------------------------------------------------------
@@ -76,17 +76,20 @@
 void
 pylith::faults::TestLiuCosSlipFn::testDbFinalSlip(void)
 { // testDbFinalSlip
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABC";
   LiuCosSlipFn slipfn;
   
   spatialdata::spatialdb::SimpleDB db(label);
   slipfn.dbFinalSlip(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbFinalSlip);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbFinalSlip->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbSlipTime);
-  CPPUNIT_ASSERT(0 == slipfn._dbRiseTime);
+  CPPUNIT_ASSERT(slipfn._dbFinalSlip);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbFinalSlip->label()));
+  CPPUNIT_ASSERT(!slipfn._dbSlipTime);
+  CPPUNIT_ASSERT(!slipfn._dbRiseTime);
+
+  PYLITH_METHOD_END;
 } // testDbFinalSlip
 
 // ----------------------------------------------------------------------
@@ -94,17 +97,20 @@
 void
 pylith::faults::TestLiuCosSlipFn::testDbSlipTime(void)
 { // testDbSlipTime
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABCD";
   LiuCosSlipFn slipfn;
   
   spatialdata::spatialdb::SimpleDB db(label);
   slipfn.dbSlipTime(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbSlipTime);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbSlipTime->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbFinalSlip);
-  CPPUNIT_ASSERT(0 == slipfn._dbRiseTime);
+  CPPUNIT_ASSERT(slipfn._dbSlipTime);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbSlipTime->label()));
+  CPPUNIT_ASSERT(!slipfn._dbFinalSlip);
+  CPPUNIT_ASSERT(!slipfn._dbRiseTime);
+
+  PYLITH_METHOD_END;
 } // testDbSlipTime
 
 // ----------------------------------------------------------------------
@@ -112,17 +118,20 @@
 void
 pylith::faults::TestLiuCosSlipFn::testDbRiseTime(void)
 { // testDbRiseTime
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABCDE";
   LiuCosSlipFn slipfn;
   
   spatialdata::spatialdb::SimpleDB db(label);
   slipfn.dbRiseTime(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbRiseTime);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbRiseTime->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbFinalSlip);
-  CPPUNIT_ASSERT(0 == slipfn._dbSlipTime);
+  CPPUNIT_ASSERT(slipfn._dbRiseTime);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbRiseTime->label()));
+  CPPUNIT_ASSERT(!slipfn._dbFinalSlip);
+  CPPUNIT_ASSERT(!slipfn._dbSlipTime);
+
+  PYLITH_METHOD_END;
 } // testDbRiseTime
 
 // ----------------------------------------------------------------------
@@ -130,6 +139,8 @@
 void
 pylith::faults::TestLiuCosSlipFn::testInitialize1D(void)
 { // testInitialize1D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/line2.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -154,6 +165,8 @@
 				       riseTimeE,
 				       numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize1D
 
 // ----------------------------------------------------------------------
@@ -161,6 +174,8 @@
 void
 pylith::faults::TestLiuCosSlipFn::testInitialize2D(void)
 { // testInitialize2D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tri3.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -186,6 +201,8 @@
 				       riseTimeE,
 				       numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize2D
 
 // ----------------------------------------------------------------------
@@ -193,6 +210,8 @@
 void
 pylith::faults::TestLiuCosSlipFn::testInitialize3D(void)
 { // testInitialize3D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tet4.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -219,6 +238,8 @@
 				       riseTimeE,
 				       numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize3D
 
 // ----------------------------------------------------------------------
@@ -237,16 +258,9 @@
   LiuCosSlipFn slipfn;
   _initialize(&mesh, &faultMesh, &slipfn, originTime);
   
-  const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
-  CPPUNIT_ASSERT(0 != cs);
-
+  const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();CPPUNIT_ASSERT(cs);
   const int spaceDim = cs->spaceDim();
-  DM dmMesh = faultMesh.dmMesh();
-  PetscErrorCode err;
-  CPPUNIT_ASSERT(dmMesh);
 
-  PetscInt       vStart, vEnd;
-  err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
   topology::Field<topology::SubMesh> slip(faultMesh);
   slip.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
   slip.allocate();
@@ -254,33 +268,35 @@
   const PylithScalar t = 2.134;
   slipfn.slip(&slip, originTime+t);
 
+  PetscDM dmMesh = faultMesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
+
+  topology::VecVisitorMesh slipVisitor(slip);
+  const PetscScalar* slipArray = slipVisitor.localArray();CPPUNIT_ASSERT(slipArray);
+
   const PylithScalar tolerance = 1.0e-06;
-  int iPoint = 0;
-  PetscSection slipSection = slip.petscSection();
-  Vec          slipVec     = slip.localVector();
-  PetscScalar *slipArray;
-  CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
-  err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
+  for(PetscInt v = vStart, iPoint = 0; v < vEnd; ++v, ++iPoint) {
     PylithScalar slipMag = 0.0;
-    for (int iDim=0; iDim < spaceDim; ++iDim)
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
       slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
+    } // for
     slipMag = sqrt(slipMag);
 
     const PylithScalar slipNorm = (slipMag > 0.0) ?
-      _slipFn(t - slipTimeE[iPoint], slipMag, riseTimeE[iPoint]) / slipMag :
-      0.0;
-    PetscInt dof, off;
+      _slipFn(t - slipTimeE[iPoint], slipMag, riseTimeE[iPoint]) / slipMag : 0.0;
 
-    err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
+    const PetscInt off = slipVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, slipVisitor.sectionDof(v));
 
-    for(PetscInt d = 0; d < dof; ++d) {
+    for(PetscInt d = 0; d < spaceDim; ++d) {
       const PylithScalar slipE = finalSlipE[iPoint*spaceDim+d] * slipNorm;
       CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, slipArray[off+d], tolerance);
     } // for
   } // for
+
+  PYLITH_METHOD_END;
 } // testSlip
 
 // ----------------------------------------------------------------------
@@ -288,6 +304,8 @@
 void
 pylith::faults::TestLiuCosSlipFn::testSlipTH(void)
 { // testSlipTH
+  PYLITH_METHOD_BEGIN;
+
   const PylithScalar t = 0.734;
   const PylithScalar finalSlip = 4.64;
   const PylithScalar riseTime = 3.23;
@@ -304,6 +322,8 @@
 
   slip = LiuCosSlipFn::_slipFn(1.0e+10, finalSlip, riseTime);
   CPPUNIT_ASSERT_DOUBLES_EQUAL(finalSlip, slip, tolerance);
+
+  PYLITH_METHOD_END;
 } // testSlipTH
 
 // ----------------------------------------------------------------------
@@ -314,7 +334,9 @@
 					      LiuCosSlipFn* slipfn,
 					      const PylithScalar originTime)
 { // _initialize
-  assert(0 != slipfn);
+  PYLITH_METHOD_BEGIN;
+
+  assert(slipfn);
   PetscErrorCode  err;
 
   const char* meshFilename = "data/tri3.mesh";
@@ -338,32 +360,8 @@
   mesh->coordsys(&cs);
 
   // Create fault mesh
-  DM       dmMesh = mesh->dmMesh(), faultBoundaryDM;
-  PetscInt firstFaultVertex = 0;
-  PetscInt firstLagrangeVertex, firstFaultCell;
-  DMLabel  groupField;
-  const bool useLagrangeConstraints = true;
+  TestSlipFn::_createFaultMesh(faultMesh, mesh, faultLabel, faultId);
 
-  err = DMPlexGetStratumSize(dmMesh, faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
-  firstFaultCell = firstLagrangeVertex;
-  if (useLagrangeConstraints) {
-    firstFaultCell += firstLagrangeVertex;
-  }
-  err = DMPlexGetLabel(dmMesh, faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  CohesiveTopology::createFault(faultMesh, faultBoundary, faultBoundaryDM,
-                                *mesh, groupField);
-  CohesiveTopology::create(mesh, *faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           faultId, firstFaultVertex, firstLagrangeVertex,
-			   firstFaultCell, useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are not
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(mesh, faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
   spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -388,6 +386,8 @@
   slipfn->dbRiseTime(&dbRiseTime);
   
   slipfn->initialize(*faultMesh, normalizer, originTime);
+
+  PYLITH_METHOD_END;
 } // _initialize
 
 // ----------------------------------------------------------------------
@@ -395,6 +395,8 @@
 void
 pylith::faults::TestLiuCosSlipFn::_testInitialize(const _TestLiuCosSlipFn::DataStruct& data)
 { // _testInitialize
+  PYLITH_METHOD_BEGIN;
+
   PetscErrorCode err;
   // Setup mesh
   topology::Mesh mesh;
@@ -413,32 +415,8 @@
 
   // Create fault mesh
   topology::SubMesh faultMesh;
-  DM       dmMesh = mesh.dmMesh(), faultBoundaryDM;
-  PetscInt firstFaultVertex = 0;
-  PetscInt firstLagrangeVertex, firstFaultCell;
-  DMLabel  groupField;
-  const bool useLagrangeConstraints = true;
+  TestSlipFn::_createFaultMesh(&faultMesh, &mesh, data.faultLabel, data.faultId);
 
-  err = DMPlexGetStratumSize(dmMesh, data.faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
-  firstFaultCell = firstLagrangeVertex;
-  if (useLagrangeConstraints) {
-    firstFaultCell += firstLagrangeVertex;
-  }
-  err = DMPlexGetLabel(dmMesh, data.faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  CohesiveTopology::createFault(&faultMesh, faultBoundary, faultBoundaryDM,
-                                mesh, groupField);
-  CohesiveTopology::create(&mesh, faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           data.faultId, firstFaultVertex, firstLagrangeVertex,
-			   firstFaultCell, useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are not
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(&mesh, &faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
   spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -466,48 +444,40 @@
   
   slipfn.initialize(faultMesh, normalizer, originTime);
 
-  CPPUNIT_ASSERT(0 != slipfn._parameters);
-  PetscSection finalSlipSection = slipfn._parameters->get("final slip").petscSection();
-  Vec          finalSlipVec     = slipfn._parameters->get("final slip").localVector();
-  PetscScalar *finalSlipArray;
-  CPPUNIT_ASSERT(finalSlipSection);CPPUNIT_ASSERT(finalSlipVec);
-  PetscSection slipTimeSection = slipfn._parameters->get("slip time").petscSection();
-  Vec          slipTimeVec     = slipfn._parameters->get("slip time").localVector();
-  PetscScalar *slipTimeArray;
-  CPPUNIT_ASSERT(slipTimeSection);CPPUNIT_ASSERT(slipTimeVec);
-  PetscSection riseTimeSection = slipfn._parameters->get("rise time").petscSection();
-  Vec          riseTimeVec     = slipfn._parameters->get("rise time").localVector();
-  PetscScalar *riseTimeArray;
-  CPPUNIT_ASSERT(riseTimeSection);CPPUNIT_ASSERT(riseTimeVec);
+  CPPUNIT_ASSERT(slipfn._parameters);
+  topology::VecVisitorMesh finalSlipVisitor(slipfn._parameters->get("final slip"));
+  const PetscScalar* finalSlipArray = finalSlipVisitor.localArray();CPPUNIT_ASSERT(finalSlipArray);
 
+  topology::VecVisitorMesh slipTimeVisitor(slipfn._parameters->get("slip time"));
+  const PetscScalar* slipTimeArray = slipTimeVisitor.localArray();CPPUNIT_ASSERT(slipTimeArray);
+
+  topology::VecVisitorMesh riseTimeVisitor(slipfn._parameters->get("rise time"));
+  const PetscScalar* riseTimeArray = riseTimeVisitor.localArray();CPPUNIT_ASSERT(riseTimeArray);
+
+  PetscDM faultDMMesh = faultMesh.dmMesh();CPPUNIT_ASSERT(faultDMMesh);
+  topology::Stratum verticesStratum(faultDMMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
+
   const PylithScalar tolerance = 1.0e-06;
-  PetscInt           vStart, vEnd, iPoint = 0;
-  err = DMPlexGetDepthStratum(faultMesh.dmMesh(), 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(slipTimeVec,  &slipTimeArray);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(riseTimeVec,  &riseTimeArray);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
-    PetscInt fsdof, fsoff, stdof, stoff, rtdof, rtoff;
+  for(PetscInt v = vStart, iPoint = 0; v < vEnd; ++v, ++iPoint) {
+    const PetscInt fsoff = finalSlipVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, finalSlipVisitor.sectionDof(v));
 
-    err = PetscSectionGetDof(finalSlipSection, v, &fsdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(finalSlipSection, v, &fsoff);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetDof(riseTimeSection, v, &rtdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(riseTimeSection, v, &rtoff);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(spaceDim, fsdof);
-    for(PetscInt d = 0; d < fsdof; ++d)
-      CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+d], finalSlipArray[fsoff+d], tolerance);
+    const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(1, slipTimeVisitor.sectionDof(v));
 
-    CPPUNIT_ASSERT_EQUAL(1, stdof);
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
+    const PetscInt rtoff = riseTimeVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(1, riseTimeVisitor.sectionDof(v));
 
-    CPPUNIT_ASSERT_EQUAL(1, rtdof);
+    for(PetscInt d = 0; d < spaceDim; ++d) {
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+d], finalSlipArray[fsoff+d], tolerance);
+    } // for
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
     CPPUNIT_ASSERT_DOUBLES_EQUAL(data.riseTimeE[iPoint], riseTimeArray[rtoff], tolerance);
   } // for
-  err = VecRestoreArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(slipTimeVec,  &slipTimeArray);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(riseTimeVec,  &riseTimeArray);CHECK_PETSC_ERROR(err);
+
+  PYLITH_METHOD_END;
 } // _testInitialize
 
 // ----------------------------------------------------------------------
@@ -517,6 +487,8 @@
 					  const PylithScalar finalSlip,
 					  const PylithScalar riseTime)
 { // _slipFn
+  PYLITH_METHOD_BEGIN;
+
   const float tau = riseTime * 1.525;
   const float tau1 = 0.13 * tau;
   const float tau2 = tau - tau1;
@@ -541,69 +513,8 @@
     slip = 1.0;
   slip *= finalSlip;
 
-  return slip;
+  PYLITH_METHOD_RETURN(slip);
 } // _slipFn
 
-// ----------------------------------------------------------------------
-// Setup fault coordinates
-void
-pylith::faults::TestLiuCosSlipFn::_setupFaultCoordinates(topology::Mesh *mesh, topology::SubMesh *faultMesh)
-{ // _setupFaultCoordinates
-  const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
-  if (!faultSieveMesh.isNull()) {
-    const ALE::Obj<RealSection>& oldCoordSection = mesh->sieveMesh()->getRealSection("coordinates");
-    faultSieveMesh->setRealSection("coordinates", oldCoordSection);
-  }
 
-  DM              dmMesh      = mesh->dmMesh();
-  DM              faultDMMesh = faultMesh->dmMesh();
-  const PetscInt  spaceDim    = mesh->dimension();
-  IS              subpointIS;
-  const PetscInt *points;
-  PetscSection    coordSection, fcoordSection;
-  PetscInt        vStart, vEnd, ffStart, ffEnd;
-  PetscErrorCode  err;
-
-  CPPUNIT_ASSERT(dmMesh);
-  CPPUNIT_ASSERT(faultDMMesh);
-  err = DMPlexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetHeightStratum(faultDMMesh, 1, &ffStart, &ffEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexCreateSubpointIS(faultDMMesh, &subpointIS);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(faultDMMesh, &fcoordSection);CHECK_PETSC_ERROR(err);
-  err = PetscSectionSetChart(fcoordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    err = PetscSectionSetDof(fcoordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
-  }
-  err = PetscSectionSetUp(fcoordSection);CHECK_PETSC_ERROR(err);
-  Vec          coordVec, fcoordVec;
-  PetscScalar *coords,  *fcoords;
-  PetscInt     coordSize;
-
-  err = PetscSectionGetStorageSize(fcoordSection, &coordSize);CHECK_PETSC_ERROR(err);
-  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
-  err = VecCreate(mesh->comm(), &fcoordVec);CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(fcoordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(fcoordVec);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    PetscInt off, foff;
-
-    // Notice that subpointMap[] does not account for cohesive cells
-    err = PetscSectionGetOffset(coordSection, points[v]+(ffEnd-ffStart), &off);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(fcoordSection, v, &foff);CHECK_PETSC_ERROR(err);
-    for(PetscInt d = 0; d < spaceDim; ++d) {
-      fcoords[foff+d] = coords[off+d];
-    }
-  }
-  err = ISRestoreIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  err = DMSetCoordinatesLocal(faultDMMesh, fcoordVec);CHECK_PETSC_ERROR(err);
-} // _setupFaultCoordinates
-
-
 // End of file 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.hh	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.hh	2013-04-04 18:29:46 UTC (rev 21720)
@@ -27,11 +27,11 @@
 #if !defined(pylith_faults_testliucosslipfn_hh)
 #define pylith_faults_testliucosslipfn_hh
 
+#include "TestSlipFn.hh" // ISA TestSlipFn
+
 #include "pylith/faults/faultsfwd.hh" // USES LiuCosSlipFn
 #include "pylith/topology/topologyfwd.hh" // USES Mesh, SubMesh
 
-#include <cppunit/extensions/HelperMacros.h>
-
 /// Namespace for pylith package
 namespace pylith {
   namespace faults {
@@ -44,7 +44,7 @@
 } // pylith
 
 /// C++ unit testing for LiuCosSlipFn
-class pylith::faults::TestLiuCosSlipFn : public CppUnit::TestFixture
+class pylith::faults::TestLiuCosSlipFn : public TestSlipFn
 { // class TestLiuCosSlipFn
 
   // CPPUNIT TEST SUITE /////////////////////////////////////////////////

Added: short/3D/PyLith/trunk/unittests/libtests/faults/TestSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestSlipFn.cc	                        (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestSlipFn.cc	2013-04-04 18:29:46 UTC (rev 21720)
@@ -0,0 +1,118 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestSlipFn.hh" // Implementation of class methods
+
+#include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::TestSlipFn::_createFaultMesh(topology::SubMesh* faultMesh,
+					     topology::Mesh* mesh,
+					     const char* faultLabel,
+					     const int faultId)
+{ // _createFaultMesh
+  PYLITH_METHOD_BEGIN;
+  
+  CPPUNIT_ASSERT(faultMesh);
+  CPPUNIT_ASSERT(mesh);
+
+  PetscErrorCode err = 0;
+  { // Create mesh
+    PetscInt firstFaultVertex = 0;
+    PetscInt firstLagrangeVertex = 0, firstFaultCell = 0;
+    PetscDMLabel groupField = NULL;
+    const bool useLagrangeConstraints = true;
+    PetscDM faultBoundaryDM = NULL;
+    PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
+    
+    err = DMPlexGetStratumSize(dmMesh, faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
+    firstFaultCell = firstLagrangeVertex;
+    if (useLagrangeConstraints) {
+      firstFaultCell += firstLagrangeVertex;
+    } // if
+    err = DMPlexGetLabel(dmMesh, faultLabel, &groupField);CHECK_PETSC_ERROR(err);
+    CPPUNIT_ASSERT(groupField);
+    ALE::Obj<SieveFlexMesh> faultBoundary = 0;
+    const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+    CPPUNIT_ASSERT(!sieveMesh.isNull());
+    CohesiveTopology::createFault(faultMesh, faultBoundary, faultBoundaryDM, *mesh, groupField);
+    CohesiveTopology::create(mesh, *faultMesh, faultBoundary, faultBoundaryDM, groupField, faultId, firstFaultVertex, firstLagrangeVertex, firstFaultCell, useLagrangeConstraints);
+  } // Create mesh
+
+  { // Need to copy coordinates from mesh to fault mesh since we are not
+    // using create() instead of createParallel().
+    PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
+    PetscDM faultDMMesh = faultMesh->dmMesh();CPPUNIT_ASSERT(faultDMMesh);
+    const int  spaceDim = mesh->dimension();
+    PetscIS subpointIS = NULL;
+    const PetscInt *points = NULL;
+    PetscSection coordSection = NULL, fcoordSection = NULL;
+    PetscInt vStart, vEnd, ffStart, ffEnd;
+    
+    err = DMPlexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+    err = DMPlexGetHeightStratum(faultDMMesh, 1, &ffStart, &ffEnd);CHECK_PETSC_ERROR(err);
+    err = DMPlexCreateSubpointIS(faultDMMesh, &subpointIS);CHECK_PETSC_ERROR(err);
+    err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+    err = DMPlexGetCoordinateSection(faultDMMesh, &fcoordSection);CHECK_PETSC_ERROR(err);
+    err = PetscSectionSetChart(fcoordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
+    for(PetscInt v = vStart; v < vEnd; ++v) {
+      err = PetscSectionSetDof(fcoordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+    } // for
+    err = PetscSectionSetUp(fcoordSection);CHECK_PETSC_ERROR(err);
+    PetscVec coordVec, fcoordVec;
+    PetscScalar *coords,  *fcoords;
+    PetscInt coordSize;
+    
+    err = PetscSectionGetStorageSize(fcoordSection, &coordSize);CHECK_PETSC_ERROR(err);
+    err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+    err = VecCreate(mesh->comm(), &fcoordVec);CHECK_PETSC_ERROR(err);
+    err = VecSetSizes(fcoordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+    err = VecSetFromOptions(fcoordVec);CHECK_PETSC_ERROR(err);
+    err = ISGetIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
+    err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+    err = VecGetArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
+    for(PetscInt v = vStart; v < vEnd; ++v) {
+      PetscInt off, foff;
+      
+      // Notice that subpointMap[] does not account for cohesive cells
+      err = PetscSectionGetOffset(coordSection, points[v]+(ffEnd-ffStart), &off);CHECK_PETSC_ERROR(err);
+      err = PetscSectionGetOffset(fcoordSection, v, &foff);CHECK_PETSC_ERROR(err);
+      for(PetscInt d = 0; d < spaceDim; ++d) {
+	fcoords[foff+d] = coords[off+d];
+      } // for
+    } // for
+    err = ISRestoreIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
+    err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
+    err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+    err = VecRestoreArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
+    err = DMSetCoordinatesLocal(faultDMMesh, fcoordVec);CHECK_PETSC_ERROR(err);
+  } // Copy coordiantes
+
+
+  PYLITH_METHOD_END;
+} // _createFaultMesh

Added: short/3D/PyLith/trunk/unittests/libtests/faults/TestSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestSlipFn.hh	                        (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestSlipFn.hh	2013-04-04 18:29:46 UTC (rev 21720)
@@ -0,0 +1,134 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/faults/TestSlipFn.hh
+ *
+ * @brief C++ TestBruneSlipFn object
+ *
+ * C++ unit testing for SlipFn.
+ */
+
+#if !defined(pylith_faults_testslipfn_hh)
+#define pylith_faults_testslipfn_hh
+
+#include "pylith/topology/topologyfwd.hh" // USES Mesh, SubMesh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/// Namespace for pylith package
+namespace pylith {
+  namespace faults {
+    class TestSlipFn;
+  } // faults
+} // pylith
+
+/// C++ unit testing for SlipFn
+class pylith::faults::TestSlipFn : public CppUnit::TestFixture
+{ // class TestSlipFn
+
+  // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+  /** Adjust topology of domain mesh and create fault mesh.
+   *
+   * @param faultMesh Fault mesh.
+   * @param mesh Domain mesh.
+   * @param faultLabel Label for fault.
+   * @param faultId Material id for fault.
+   */
+  static
+  void
+  _createFaultMesh(topology::SubMesh* faultMesh,
+		   topology::Mesh* mesh,
+		   const char* faultLabel,
+		   const int faultId);
+
+}; // class TestSlipFn
+
+#endif // pylith_faults_testslipfn_hh
+
+
+// End of file 
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/faults/TestSlipFn.hh
+ *
+ * @brief C++ TestBruneSlipFn object
+ *
+ * C++ unit testing for SlipFn.
+ */
+
+#if !defined(pylith_faults_testslipfn_hh)
+#define pylith_faults_testslipfn_hh
+
+#include "pylith/topology/topologyfwd.hh" // USES Mesh, SubMesh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/// Namespace for pylith package
+namespace pylith {
+  namespace faults {
+    class TestSlipFn;
+  } // faults
+} // pylith
+
+/// C++ unit testing for SlipFn
+class pylith::faults::TestSlipFn : public CppUnit::TestFixture
+{ // class TestSlipFn
+
+  // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+  /** Adjust topology of domain mesh and create fault mesh.
+   *
+   * @param faultMesh Fault mesh.
+   * @param mesh Domain mesh.
+   * @param faultLabel Label for fault.
+   * @param faultId Material id for fault.
+   */
+  static
+  void
+  _createFaultMesh(topology::SubMesh* faultMesh,
+		   topology::Mesh* mesh,
+		   const char* faultLabel,
+		   const int faultId);
+
+}; // class TestSlipFn
+
+#endif // pylith_faults_testslipfn_hh
+
+
+// End of file 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc	2013-04-04 18:29:46 UTC (rev 21720)
@@ -22,11 +22,12 @@
 
 #include "pylith/faults/StepSlipFn.hh" // USES StepSlipFn
 
-#include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 
 #include "spatialdata/geocoords/CSCart.hh" // USES CSCart
@@ -57,16 +58,15 @@
 } // pylith
 
 // ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
 // Test constructor.
 void
 pylith::faults::TestStepSlipFn::testConstructor(void)
 { // testConstructor
+  PYLITH_METHOD_BEGIN;
+
   StepSlipFn slipfn;
+
+  PYLITH_METHOD_END;
 } // testConstructor
 
 // ----------------------------------------------------------------------
@@ -74,16 +74,19 @@
 void
 pylith::faults::TestStepSlipFn::testDbFinalSlip(void)
 { // testDbFinalSlip
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABC";
   StepSlipFn slipfn;
   
   spatialdata::spatialdb::SimpleDB db(label);
   slipfn.dbFinalSlip(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbFinalSlip);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbFinalSlip->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbSlipTime);
+  CPPUNIT_ASSERT(slipfn._dbFinalSlip);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbFinalSlip->label()));
+  CPPUNIT_ASSERT(!slipfn._dbSlipTime);
+
+  PYLITH_METHOD_END;
 } // testDbFinalSlip
 
 // ----------------------------------------------------------------------
@@ -91,16 +94,19 @@
 void
 pylith::faults::TestStepSlipFn::testDbSlipTime(void)
 { // testDbSlipTime
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABCD";
   StepSlipFn slipfn;
   
   spatialdata::spatialdb::SimpleDB db(label);
   slipfn.dbSlipTime(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbSlipTime);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbSlipTime->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbFinalSlip);
+  CPPUNIT_ASSERT(slipfn._dbSlipTime);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbSlipTime->label()));
+  CPPUNIT_ASSERT(!slipfn._dbFinalSlip);
+
+  PYLITH_METHOD_END;
 } // testDbSlipTime
 
 // ----------------------------------------------------------------------
@@ -108,6 +114,8 @@
 void
 pylith::faults::TestStepSlipFn::testInitialize1D(void)
 { // testInitialize1D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/line2.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -128,6 +136,8 @@
 				      slipTimeE,
 				      numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize1D
 
 // ----------------------------------------------------------------------
@@ -135,6 +145,8 @@
 void
 pylith::faults::TestStepSlipFn::testInitialize2D(void)
 { // testInitialize2D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tri3.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -156,6 +168,8 @@
 				      slipTimeE,
 				      numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize2D
 
 // ----------------------------------------------------------------------
@@ -163,6 +177,8 @@
 void
 pylith::faults::TestStepSlipFn::testInitialize3D(void)
 { // testInitialize3D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tet4.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -185,6 +201,8 @@
 				      slipTimeE,
 				      numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize3D
 
 // ----------------------------------------------------------------------
@@ -192,6 +210,8 @@
 void
 pylith::faults::TestStepSlipFn::testSlip(void)
 { // testSlip
+  PYLITH_METHOD_BEGIN;
+
   const PylithScalar slipE[] = { 2.3, 0.1, 
 			   0.0, 0.0};
   const PylithScalar originTime = 5.064;
@@ -201,16 +221,9 @@
   StepSlipFn slipfn;
   _initialize(&mesh, &faultMesh, &slipfn, originTime);
   
-  const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
-  CPPUNIT_ASSERT(0 != cs);
-
+  const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();CPPUNIT_ASSERT(cs);
   const int spaceDim = cs->spaceDim();
-  DM dmMesh = faultMesh.dmMesh();
-  PetscErrorCode err;
-  CPPUNIT_ASSERT(dmMesh);
 
-  PetscInt       vStart, vEnd;
-  err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
   topology::Field<topology::SubMesh> slip(faultMesh);
   slip.newSection(topology::Field<topology::Mesh>::VERTICES_FIELD, spaceDim);
   slip.allocate();
@@ -218,23 +231,25 @@
   const PylithScalar t = 1.234;
   slipfn.slip(&slip, originTime+t);
 
+  PetscDM dmMesh = faultMesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
+
+  topology::VecVisitorMesh slipVisitor(slip);
+  const PetscScalar* slipArray = slipVisitor.localArray();CPPUNIT_ASSERT(slipArray);
+
   const PylithScalar tolerance = 1.0e-06;
-  PetscScalar       *slipArray;
-  int iPoint = 0;
-  PetscSection slipSection = slip.petscSection();
-  Vec          slipVec     = slip.localVector();
-  CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
-  err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
-    PetscInt dof, off;
+  for(PetscInt v = vStart, iPoint = 0; v < vEnd; ++v, ++iPoint) {
 
-    err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
-    for(PetscInt d = 0; d < dof; ++d)
+    const PetscInt off = slipVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, slipVisitor.sectionDof(v));
+
+    for(PetscInt d = 0; d < spaceDim; ++d)
       CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE[iPoint*spaceDim+d], slipArray[off+d], tolerance);
   } // for
-  err = VecRestoreArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+
+  PYLITH_METHOD_END;
 } // testSlip
 
 // ----------------------------------------------------------------------
@@ -245,9 +260,11 @@
 					    StepSlipFn* slipfn,
 					    const PylithScalar originTime)
 { // _initialize
-  CPPUNIT_ASSERT(0 != mesh);
-  CPPUNIT_ASSERT(0 != faultMesh);
-  CPPUNIT_ASSERT(0 != slipfn);
+  PYLITH_METHOD_BEGIN;
+
+  CPPUNIT_ASSERT(mesh);
+  CPPUNIT_ASSERT(faultMesh);
+  CPPUNIT_ASSERT(slipfn);
   PetscErrorCode err;
 
   const char* meshFilename = "data/tri3.mesh";
@@ -270,33 +287,8 @@
   mesh->coordsys(&cs);
 
   // Create fault mesh
-  DM       dmMesh = mesh->dmMesh(), faultBoundaryDM;
-  PetscInt firstFaultVertex = 0;
-  PetscInt firstLagrangeVertex, firstFaultCell;
-  DMLabel  groupField;
-  const bool useLagrangeConstraints = true;
+  TestSlipFn::_createFaultMesh(faultMesh, mesh, faultLabel, faultId);
 
-  err = DMPlexGetStratumSize(dmMesh, faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
-  firstFaultCell = firstLagrangeVertex;
-  if (useLagrangeConstraints) {
-    firstFaultCell += firstLagrangeVertex;
-  }
-  err = DMPlexGetLabel(dmMesh, faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());CPPUNIT_ASSERT(dmMesh);
-  CohesiveTopology::createFault(faultMesh, faultBoundary, faultBoundaryDM,
-                                *mesh, groupField);
-  CohesiveTopology::create(mesh, *faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           faultId,
-                           firstFaultVertex, firstLagrangeVertex, firstFaultCell,
-                           useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are not
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(mesh, faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
   spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -315,6 +307,8 @@
   slipfn->dbSlipTime(&dbSlipTime);
   
   slipfn->initialize(*faultMesh, normalizer, originTime);
+
+  PYLITH_METHOD_END;
 } // _initialize
 
 // ----------------------------------------------------------------------
@@ -322,6 +316,8 @@
 void
 pylith::faults::TestStepSlipFn::_testInitialize(const _TestStepSlipFn::DataStruct& data)
 { // _testInitialize
+  PYLITH_METHOD_BEGIN;
+
   PetscErrorCode err;
   // Setup mesh
   topology::Mesh mesh;
@@ -340,32 +336,8 @@
 
   // Create fault mesh
   topology::SubMesh faultMesh;
-  int firstFaultVertex    = 0;
-  int firstLagrangeVertex = mesh.sieveMesh()->getIntSection(data.faultLabel)->size();
-  int firstFaultCell      = mesh.sieveMesh()->getIntSection(data.faultLabel)->size();
-  const bool useLagrangeConstraints = true;
-  if (useLagrangeConstraints) {
-    firstFaultCell += mesh.sieveMesh()->getIntSection(data.faultLabel)->size();
-  }
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  DM      dmMesh = mesh.dmMesh(), faultBoundaryDM;
-  DMLabel groupField;
+  TestSlipFn::_createFaultMesh(&faultMesh, &mesh, data.faultLabel, data.faultId);
 
-  CPPUNIT_ASSERT(!sieveMesh.isNull());CPPUNIT_ASSERT(dmMesh);
-  err = DMPlexGetLabel(dmMesh, data.faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  CohesiveTopology::createFault(&faultMesh, faultBoundary, faultBoundaryDM,
-                                mesh, groupField);
-  CohesiveTopology::create(&mesh, faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           data.faultId,
-                           firstFaultVertex, firstLagrangeVertex, firstFaultCell,
-                           useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(&mesh, &faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
   spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -387,105 +359,34 @@
   
   slipfn.initialize(faultMesh, normalizer, originTime);
 
-  CPPUNIT_ASSERT(0 != slipfn._parameters);
-  PetscSection finalSlipSection = slipfn._parameters->get("final slip").petscSection();
-  Vec          finalSlipVec     = slipfn._parameters->get("final slip").localVector();
-  CPPUNIT_ASSERT(finalSlipSection);CPPUNIT_ASSERT(finalSlipVec);
-  PetscSection slipTimeSection = slipfn._parameters->get("slip time").petscSection();
-  Vec          slipTimeVec     = slipfn._parameters->get("slip time").localVector();
-  CPPUNIT_ASSERT(slipTimeSection);CPPUNIT_ASSERT(slipTimeVec);
+  CPPUNIT_ASSERT(slipfn._parameters);
+  topology::VecVisitorMesh finalSlipVisitor(slipfn._parameters->get("final slip"));
+  const PetscScalar* finalSlipArray = finalSlipVisitor.localArray();CPPUNIT_ASSERT(finalSlipArray);
 
-  const PylithScalar tolerance = 1.0e-06;
-  PetscScalar       *finalSlipArray, *slipTimeArray;
-  PetscInt           vStart, vEnd, iPoint = 0;
-  err = DMPlexGetDepthStratum(faultMesh.dmMesh(), 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
-    PetscInt dof, off;
+  topology::VecVisitorMesh slipTimeVisitor(slipfn._parameters->get("slip time"));
+  const PetscScalar* slipTimeArray = slipTimeVisitor.localArray();CPPUNIT_ASSERT(slipTimeArray);
 
-    err = PetscSectionGetDof(finalSlipSection, v, &dof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(finalSlipSection, v, &off);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
-    for(PetscInt d = 0; d < spaceDim; ++d)
-      CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+d], finalSlipArray[off+d], tolerance);
+  PetscDM faultDMMesh = faultMesh.dmMesh();CPPUNIT_ASSERT(faultDMMesh);
+  topology::Stratum verticesStratum(faultDMMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
 
-    err = PetscSectionGetDof(slipTimeSection, v, &dof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipTimeSection, v, &off);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(1, dof);
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[off], tolerance);
-  } // for
-  err = VecRestoreArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
-} // _testInitialize
+  const PylithScalar tolerance = 1.0e-06;
+  for(PetscInt v = vStart, iPoint = 0; v < vEnd; ++v, ++iPoint) {
+    const PetscInt fsoff = finalSlipVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, finalSlipVisitor.sectionDof(v));
 
-// ----------------------------------------------------------------------
-// Setup fault coordinates
-void
-pylith::faults::TestStepSlipFn::_setupFaultCoordinates(topology::Mesh *mesh, topology::SubMesh *faultMesh)
-{ // _setupFaultCoordinates
-  const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
-  if (!faultSieveMesh.isNull()) {
-    const ALE::Obj<RealSection>& oldCoordSection = mesh->sieveMesh()->getRealSection("coordinates");
-    faultSieveMesh->setRealSection("coordinates", oldCoordSection);
-  }
+    const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(1, slipTimeVisitor.sectionDof(v));
 
-  DM              dmMesh      = mesh->dmMesh();
-  DM              faultDMMesh = faultMesh->dmMesh();
-  const PetscInt  spaceDim    = mesh->dimension();
-  IS              subpointIS;
-  const PetscInt *points;
-  PetscSection    coordSection, fcoordSection;
-  PetscInt        vStart, vEnd, ffStart, ffEnd;
-  PetscErrorCode  err;
-
-  CPPUNIT_ASSERT(dmMesh);
-  CPPUNIT_ASSERT(faultDMMesh);
-  err = DMPlexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetHeightStratum(faultDMMesh, 1, &ffStart, &ffEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexCreateSubpointIS(faultDMMesh, &subpointIS);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(faultDMMesh, &fcoordSection);CHECK_PETSC_ERROR(err);
-  err = PetscSectionSetChart(fcoordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    err = PetscSectionSetDof(fcoordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
-  }
-  err = PetscSectionSetUp(fcoordSection);CHECK_PETSC_ERROR(err);
-  Vec          coordVec, fcoordVec;
-  PetscScalar *coords,  *fcoords;
-  PetscInt     coordSize;
-
-  err = PetscSectionGetStorageSize(fcoordSection, &coordSize);CHECK_PETSC_ERROR(err);
-  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
-  err = VecCreate(mesh->comm(), &fcoordVec);CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(fcoordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(fcoordVec);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    PetscInt off, foff;
-
-    // Notice that subpointMap[] does not account for cohesive cells
-    err = PetscSectionGetOffset(coordSection, points[v]+(ffEnd-ffStart), &off);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(fcoordSection, v, &foff);CHECK_PETSC_ERROR(err);
     for(PetscInt d = 0; d < spaceDim; ++d) {
-      fcoords[foff+d] = coords[off+d];
-    }
-    if (!faultSieveMesh.isNull()) {
-      const PetscScalar *oldCoords = mesh->sieveMesh()->getRealSection("coordinates")->restrictPoint(points[v]);
-      for(PetscInt d = 0; d < spaceDim; ++d) {
-        CPPUNIT_ASSERT_DOUBLES_EQUAL(oldCoords[d], coords[off+d], 1.0e-6);
-      }
-    }
-  }
-  err = ISRestoreIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  err = DMSetCoordinatesLocal(faultDMMesh, fcoordVec);CHECK_PETSC_ERROR(err);
-  err = VecDestroy(&fcoordVec);CHECK_PETSC_ERROR(err);
-} // _setupFaultCoordinates
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+d], finalSlipArray[fsoff+d], tolerance);
+    } // for
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
+  } // for
 
+  PYLITH_METHOD_END;
+} // _testInitialize
 
+
 // End of file 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.hh	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.hh	2013-04-04 18:29:46 UTC (rev 21720)
@@ -27,11 +27,11 @@
 #if !defined(pylith_faults_teststepslipfn_hh)
 #define pylith_faults_teststepslipfn_hh
 
+#include "TestSlipFn.hh" // ISA TestSlipFn
+
 #include "pylith/faults/faultsfwd.hh" // USES StepSlipFn
 #include "pylith/topology/topologyfwd.hh" // USES Mesh, SubMesh
 
-#include <cppunit/extensions/HelperMacros.h>
-
 /// Namespace for pylith package
 namespace pylith {
   namespace faults {
@@ -44,7 +44,7 @@
 } // pylith
 
 /// C++ unit testing for StepSlipFn
-class pylith::faults::TestStepSlipFn : public CppUnit::TestFixture
+class pylith::faults::TestStepSlipFn : public TestSlipFn
 { // class TestStepSlipFn
 
   // CPPUNIT TEST SUITE /////////////////////////////////////////////////

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.cc	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.cc	2013-04-04 18:29:46 UTC (rev 21720)
@@ -22,11 +22,12 @@
 
 #include "pylith/faults/TimeHistorySlipFn.hh" // USES TimeHistorySlipFn
 
-#include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 
 #include "spatialdata/geocoords/CSCart.hh" // USES CSCart
@@ -39,11 +40,6 @@
 CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestTimeHistorySlipFn );
 
 // ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
 namespace pylith {
   namespace faults {
     namespace _TestTimeHistorySlipFn {
@@ -68,7 +64,11 @@
 void
 pylith::faults::TestTimeHistorySlipFn::testConstructor(void)
 { // testConstructor
+  PYLITH_METHOD_BEGIN;
+
   TimeHistorySlipFn slipfn;
+
+  PYLITH_METHOD_END;
 } // testConstructor
 
 // ----------------------------------------------------------------------
@@ -76,17 +76,20 @@
 void
 pylith::faults::TestTimeHistorySlipFn::testDbAmplitude(void)
 { // testDbAmplitude
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABC";
   TimeHistorySlipFn slipfn;
   
   spatialdata::spatialdb::SimpleDB db(label);
   slipfn.dbAmplitude(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbAmplitude);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbAmplitude->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbSlipTime);
-  CPPUNIT_ASSERT(0 == slipfn._dbTimeHistory);
+  CPPUNIT_ASSERT(slipfn._dbAmplitude);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbAmplitude->label()));
+  CPPUNIT_ASSERT(!slipfn._dbSlipTime);
+  CPPUNIT_ASSERT(!slipfn._dbTimeHistory);
+
+  PYLITH_METHOD_END;
 } // testDbAmplitude
 
 // ----------------------------------------------------------------------
@@ -94,17 +97,20 @@
 void
 pylith::faults::TestTimeHistorySlipFn::testDbSlipTime(void)
 { // testDbSlipTime
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABCD";
   TimeHistorySlipFn slipfn;
   
   spatialdata::spatialdb::SimpleDB db(label);
   slipfn.dbSlipTime(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbSlipTime);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbSlipTime->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbAmplitude);
-  CPPUNIT_ASSERT(0 == slipfn._dbTimeHistory);
+  CPPUNIT_ASSERT(slipfn._dbSlipTime);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbSlipTime->label()));
+  CPPUNIT_ASSERT(!slipfn._dbAmplitude);
+  CPPUNIT_ASSERT(!slipfn._dbTimeHistory);
+
+  PYLITH_METHOD_END;
 } // testDbSlipTime
 
 // ----------------------------------------------------------------------
@@ -112,17 +118,20 @@
 void
 pylith::faults::TestTimeHistorySlipFn::testDbTimeHistory(void)
 { // testDbTimeHistory
+  PYLITH_METHOD_BEGIN;
+
   const char* label = "database ABCDE";
   TimeHistorySlipFn slipfn;
   
   spatialdata::spatialdb::TimeHistory db(label);
   slipfn.dbTimeHistory(&db);
 
-  CPPUNIT_ASSERT(0 != slipfn._dbTimeHistory);
-  CPPUNIT_ASSERT_EQUAL(std::string(label),
-		       std::string(slipfn._dbTimeHistory->label()));
-  CPPUNIT_ASSERT(0 == slipfn._dbAmplitude);
-  CPPUNIT_ASSERT(0 == slipfn._dbSlipTime);
+  CPPUNIT_ASSERT(slipfn._dbTimeHistory);
+  CPPUNIT_ASSERT_EQUAL(std::string(label), std::string(slipfn._dbTimeHistory->label()));
+  CPPUNIT_ASSERT(!slipfn._dbAmplitude);
+  CPPUNIT_ASSERT(!slipfn._dbSlipTime);
+
+  PYLITH_METHOD_END;
 } // testDbTimeHistory
 
 // ----------------------------------------------------------------------
@@ -130,6 +139,8 @@
 void
 pylith::faults::TestTimeHistorySlipFn::testInitialize1D(void)
 { // testInitialize1D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/line2.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -152,6 +163,8 @@
 				       slipTimeE,
 				       numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize1D
 
 // ----------------------------------------------------------------------
@@ -159,6 +172,8 @@
 void
 pylith::faults::TestTimeHistorySlipFn::testInitialize2D(void)
 { // testInitialize2D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tri3.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -182,6 +197,8 @@
 				       slipTimeE,
 				       numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize2D
 
 // ----------------------------------------------------------------------
@@ -189,6 +206,8 @@
 void
 pylith::faults::TestTimeHistorySlipFn::testInitialize3D(void)
 { // testInitialize3D
+  PYLITH_METHOD_BEGIN;
+
   const char* meshFilename = "data/tet4.mesh";
   const char* faultLabel = "fault";
   const int faultId = 2;
@@ -213,6 +232,8 @@
 				       slipTimeE,
 				       numConstraintPts};
   _testInitialize(data);
+
+  PYLITH_METHOD_END;
 } // testInitialize3D
 
 // ----------------------------------------------------------------------
@@ -220,6 +241,8 @@
 void
 pylith::faults::TestTimeHistorySlipFn::testSlip(void)
 { // testSlip
+  PYLITH_METHOD_BEGIN;
+
   const PylithScalar slipTimeE[] = { 1.2, 1.3 };
   const PylithScalar slipE[] = { 0.92, 0.04,
 			   0.84, 0.07 };
@@ -231,16 +254,9 @@
   spatialdata::spatialdb::TimeHistory th;
   _initialize(&mesh, &faultMesh, &slipfn, &th, originTime);
   
-  const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
-  CPPUNIT_ASSERT(0 != cs);
-
+  const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();CPPUNIT_ASSERT(cs);
   const int spaceDim = cs->spaceDim();
-  DM dmMesh = faultMesh.dmMesh();
-  PetscErrorCode err;
-  CPPUNIT_ASSERT(dmMesh);
 
-  PetscInt       vStart, vEnd;
-  err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
   topology::Field<topology::SubMesh> slip(faultMesh);
   slip.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
   slip.allocate();
@@ -248,35 +264,39 @@
   const PylithScalar t = 2.0;
   slipfn.slip(&slip, originTime+t);
 
-  const PylithScalar tolerance = 1.0e-06;
-  int iPoint = 0;
-  PetscSection slipSection = slip.petscSection();
-  Vec          slipVec     = slip.localVector();
-  PetscScalar *slipArray;
-  CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
-  err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
-    PetscInt dof, off;
+  PetscDM dmMesh = faultMesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
 
-    err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
+  topology::VecVisitorMesh slipVisitor(slip);
+  const PetscScalar* slipArray = slipVisitor.localArray();CPPUNIT_ASSERT(slipArray);
+
+  const PylithScalar tolerance = 1.0e-06;
+  for(PetscInt v = vStart, iPoint = 0; v < vEnd; ++v, ++iPoint) {
+    const PetscInt off = slipVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, slipVisitor.sectionDof(v));
     
-    for(PetscInt d = 0; d < dof; ++d)
+    for(PetscInt d = 0; d < spaceDim; ++d) {
       CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE[iPoint*spaceDim+d], slipArray[off+d], tolerance);
+    } // for
   } // for
+
+  PYLITH_METHOD_END;
 } // testSlip
 
 // ----------------------------------------------------------------------
 // Initialize TimeHistorySlipFn.
 void
 pylith::faults::TestTimeHistorySlipFn::_initialize(topology::Mesh* mesh,
-					      topology::SubMesh* faultMesh,
-					      TimeHistorySlipFn* slipfn,
-					      spatialdata::spatialdb::TimeHistory* th,
-					      const PylithScalar originTime)
+						   topology::SubMesh* faultMesh,
+						   TimeHistorySlipFn* slipfn,
+						   spatialdata::spatialdb::TimeHistory* th,
+						   const PylithScalar originTime)
 { // _initialize
-  assert(0 != slipfn);
+  PYLITH_METHOD_BEGIN;
+
+  assert(slipfn);
   PetscErrorCode err;
 
   const char* meshFilename = "data/tri3.mesh";
@@ -300,32 +320,8 @@
   mesh->coordsys(&cs);
 
   // Create fault mesh
-  DM       dmMesh = mesh->dmMesh(), faultBoundaryDM;
-  PetscInt firstFaultVertex = 0;
-  PetscInt firstLagrangeVertex, firstFaultCell;
-  DMLabel  groupField;
-  const bool useLagrangeConstraints = true;
+  TestSlipFn::_createFaultMesh(faultMesh, mesh, faultLabel, faultId);
 
-  err = DMPlexGetStratumSize(dmMesh, faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
-  firstFaultCell = firstLagrangeVertex;
-  if (useLagrangeConstraints) {
-    firstFaultCell += firstLagrangeVertex;
-  }
-  err = DMPlexGetLabel(dmMesh, faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  CohesiveTopology::createFault(faultMesh, faultBoundary, faultBoundaryDM,
-                                *mesh, groupField);
-  CohesiveTopology::create(mesh, *faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           faultId, firstFaultVertex, firstLagrangeVertex, 
-			   firstFaultCell, useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are not
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(mesh, faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbAmplitude("slip amplitude");
   spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -348,6 +344,8 @@
   slipfn->dbTimeHistory(th);
   
   slipfn->initialize(*faultMesh, normalizer, originTime);
+
+  PYLITH_METHOD_END;
 } // _initialize
 
 // ----------------------------------------------------------------------
@@ -355,6 +353,8 @@
 void
 pylith::faults::TestTimeHistorySlipFn::_testInitialize(const _TestTimeHistorySlipFn::DataStruct& data)
 { // _testInitialize
+  PYLITH_METHOD_BEGIN;
+
   PetscErrorCode err;
   // Setup mesh
   topology::Mesh mesh;
@@ -373,32 +373,8 @@
 
   // Create fault mesh
   topology::SubMesh faultMesh;
-  DM       dmMesh = mesh.dmMesh(), faultBoundaryDM;
-  PetscInt firstFaultVertex = 0;
-  PetscInt firstLagrangeVertex, firstFaultCell;
-  DMLabel  groupField;
-  const bool useLagrangeConstraints = true;
+  TestSlipFn::_createFaultMesh(&faultMesh, &mesh, data.faultLabel, data.faultId);
 
-  err = DMPlexGetStratumSize(dmMesh, data.faultLabel, 1, &firstLagrangeVertex);CHECK_PETSC_ERROR(err);
-  firstFaultCell = firstLagrangeVertex;
-  if (useLagrangeConstraints) {
-    firstFaultCell += firstLagrangeVertex;
-  }
-  err = DMPlexGetLabel(dmMesh, data.faultLabel, &groupField);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT(groupField);
-  ALE::Obj<SieveFlexMesh> faultBoundary = 0;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  CohesiveTopology::createFault(&faultMesh, faultBoundary, faultBoundaryDM,
-                                mesh, groupField);
-  CohesiveTopology::create(&mesh, faultMesh, faultBoundary, faultBoundaryDM,
-                           groupField,
-                           data.faultId, firstFaultVertex, firstLagrangeVertex,
-			   firstFaultCell, useLagrangeConstraints);
-  // Need to copy coordinates from mesh to fault mesh since we are not
-  // using create() instead of createParallel().
-  _setupFaultCoordinates(&mesh, &faultMesh);
-
   // Setup databases
   spatialdata::spatialdb::SimpleDB dbAmplitude("slip amplitude");
   spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -424,97 +400,34 @@
   
   slipfn.initialize(faultMesh, normalizer, originTime);
 
-  CPPUNIT_ASSERT(0 != slipfn._parameters);
-  PetscSection finalSlipSection = slipfn._parameters->get("slip amplitude").petscSection();
-  Vec          finalSlipVec     = slipfn._parameters->get("slip amplitude").localVector();
-  PetscScalar *finalSlipArray;
-  CPPUNIT_ASSERT(finalSlipSection);CPPUNIT_ASSERT(finalSlipVec);
-  PetscSection slipTimeSection = slipfn._parameters->get("slip time").petscSection();
-  Vec          slipTimeVec     = slipfn._parameters->get("slip time").localVector();
-  PetscScalar *slipTimeArray;
-  CPPUNIT_ASSERT(slipTimeSection);CPPUNIT_ASSERT(slipTimeVec);
+  CPPUNIT_ASSERT(slipfn._parameters);
+  topology::VecVisitorMesh finalSlipVisitor(slipfn._parameters->get("slip amplitude"));
+  const PetscScalar* finalSlipArray = finalSlipVisitor.localArray();CPPUNIT_ASSERT(finalSlipArray);
 
+  topology::VecVisitorMesh slipTimeVisitor(slipfn._parameters->get("slip time"));
+  const PetscScalar* slipTimeArray = slipTimeVisitor.localArray();CPPUNIT_ASSERT(slipTimeArray);
+
+  PetscDM faultDMMesh = faultMesh.dmMesh();CPPUNIT_ASSERT(faultDMMesh);
+  topology::Stratum verticesStratum(faultDMMesh, topology::Stratum::DEPTH, 0);
+  const PetscInt vStart = verticesStratum.begin();
+  const PetscInt vEnd = verticesStratum.end();
+
   const PylithScalar tolerance = 1.0e-06;
-  PetscInt           vStart, vEnd, iPoint = 0;
-  err = DMPlexGetDepthStratum(faultMesh.dmMesh(), 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(slipTimeVec,  &slipTimeArray);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
-    PetscInt fsdof, fsoff, stdof, stoff;
+  for(PetscInt v = vStart, iPoint = 0; v < vEnd; ++v, ++iPoint) {
+    const PetscInt fsoff = finalSlipVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(spaceDim, finalSlipVisitor.sectionDof(v));
 
-    err = PetscSectionGetDof(finalSlipSection, v, &fsdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(finalSlipSection, v, &fsoff);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
-    CPPUNIT_ASSERT_EQUAL(spaceDim, fsdof);
-    for(PetscInt d = 0; d < fsdof; ++d)
-      CPPUNIT_ASSERT_DOUBLES_EQUAL(data.amplitudeE[iPoint*spaceDim+d], finalSlipArray[fsoff+d], tolerance);
+    const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+    CPPUNIT_ASSERT_EQUAL(1, slipTimeVisitor.sectionDof(v));
 
-    CPPUNIT_ASSERT_EQUAL(1, stdof);
+    for(PetscInt d = 0; d < spaceDim; ++d) {
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(data.amplitudeE[iPoint*spaceDim+d], finalSlipArray[fsoff+d], tolerance);
+    } // for
     CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
   } // for
+
+  PYLITH_METHOD_END;
 } // _testInitialize
 
-// ----------------------------------------------------------------------
-// Setup fault coordinates
-void
-pylith::faults::TestTimeHistorySlipFn::_setupFaultCoordinates(topology::Mesh *mesh, topology::SubMesh *faultMesh)
-{ // _setupFaultCoordinates
-  const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
-  if (!faultSieveMesh.isNull()) {
-    const ALE::Obj<RealSection>& oldCoordSection = mesh->sieveMesh()->getRealSection("coordinates");
-    faultSieveMesh->setRealSection("coordinates", oldCoordSection);
-  }
 
-  DM              dmMesh      = mesh->dmMesh();
-  DM              faultDMMesh = faultMesh->dmMesh();
-  const PetscInt  spaceDim    = mesh->dimension();
-  IS              subpointIS;
-  const PetscInt *points;
-  PetscSection    coordSection, fcoordSection;
-  PetscInt        vStart, vEnd, ffStart, ffEnd;
-  PetscErrorCode  err;
-
-  CPPUNIT_ASSERT(dmMesh);
-  CPPUNIT_ASSERT(faultDMMesh);
-  err = DMPlexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetHeightStratum(faultDMMesh, 1, &ffStart, &ffEnd);CHECK_PETSC_ERROR(err);
-  err = DMPlexCreateSubpointIS(faultDMMesh, &subpointIS);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMPlexGetCoordinateSection(faultDMMesh, &fcoordSection);CHECK_PETSC_ERROR(err);
-  err = PetscSectionSetChart(fcoordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    err = PetscSectionSetDof(fcoordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
-  }
-  err = PetscSectionSetUp(fcoordSection);CHECK_PETSC_ERROR(err);
-  Vec          coordVec, fcoordVec;
-  PetscScalar *coords,  *fcoords;
-  PetscInt     coordSize;
-
-  err = PetscSectionGetStorageSize(fcoordSection, &coordSize);CHECK_PETSC_ERROR(err);
-  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
-  err = VecCreate(mesh->comm(), &fcoordVec);CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(fcoordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(fcoordVec);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  for(PetscInt v = vStart; v < vEnd; ++v) {
-    PetscInt off, foff;
-
-    // Notice that subpointMap[] does not account for cohesive cells
-    err = PetscSectionGetOffset(coordSection, points[v]+(ffEnd-ffStart), &off);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(fcoordSection, v, &foff);CHECK_PETSC_ERROR(err);
-    for(PetscInt d = 0; d < spaceDim; ++d) {
-      fcoords[foff+d] = coords[off+d];
-    }
-  }
-  err = ISRestoreIndices(subpointIS, &points);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
-  err = VecRestoreArray(fcoordVec, &fcoords);CHECK_PETSC_ERROR(err);
-  err = DMSetCoordinatesLocal(faultDMMesh, fcoordVec);CHECK_PETSC_ERROR(err);
-} // _setupFaultCoordinates
-
-
 // End of file 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.hh	2013-04-04 18:26:22 UTC (rev 21719)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.hh	2013-04-04 18:29:46 UTC (rev 21720)
@@ -27,12 +27,12 @@
 #if !defined(pylith_faults_testtimehistoryslipfn_hh)
 #define pylith_faults_testtimehistoryslipfn_hh
 
+#include "TestSlipFn.hh" // ISA TestSlipFn
+
 #include "pylith/faults/faultsfwd.hh" // USES TimeHistorySlipFn
 #include "pylith/topology/topologyfwd.hh" // USES Mesh, SubMesh
 #include "spatialdata/spatialdb/spatialdbfwd.hh" // USES TimeHistory
 
-#include <cppunit/extensions/HelperMacros.h>
-
 /// Namespace for pylith package
 namespace pylith {
   namespace faults {
@@ -45,7 +45,7 @@
 } // pylith
 
 /// C++ unit testing for TimeHistorySlipFn
-class pylith::faults::TestTimeHistorySlipFn : public CppUnit::TestFixture
+class pylith::faults::TestTimeHistorySlipFn : public TestSlipFn
 { // class TestTimeHistorySlipFn
 
   // CPPUNIT TEST SUITE /////////////////////////////////////////////////



More information about the CIG-COMMITS mailing list