[cig-commits] r22057 - in short/3D/PyLith/trunk: libsrc/pylith/faults libsrc/pylith/meshio libsrc/pylith/topology modulesrc/topology pylith/topology unittests/libtests/bc unittests/libtests/faults unittests/libtests/feassemble unittests/libtests/friction unittests/libtests/materials unittests/libtests/meshio unittests/libtests/topology unittests/pytests/faults unittests/pytests/topology

brad at geodynamics.org brad at geodynamics.org
Mon May 13 11:43:34 PDT 2013


Author: brad
Date: 2013-05-13 11:43:34 -0700 (Mon, 13 May 2013)
New Revision: 22057

Modified:
   short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
   short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.hh
   short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
   short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.hh
   short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc
   short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.hh
   short/3D/PyLith/trunk/libsrc/pylith/topology/ReverseCuthillMcKee.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh
   short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc
   short/3D/PyLith/trunk/modulesrc/topology/Mesh.i
   short/3D/PyLith/trunk/modulesrc/topology/MeshOps.i
   short/3D/PyLith/trunk/pylith/topology/MeshImporter.py
   short/3D/PyLith/trunk/pylith/topology/MeshImporterDist.py
   short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryConditionPoints.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBoundary.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependentPoints.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitLgDeform.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicitLgDeform.cc
   short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterFaultMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
   short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldSubMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.hh
   short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.cc
   short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.hh
   short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.hh
   short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveDyn.py
   short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py
   short/3D/PyLith/trunk/unittests/pytests/topology/TestSubMesh.py
Log:
Prepare for removal of SubMesh. Cleanup Mesh and SubMesh so we can eliminate SubMesh and templates. Refactor some Mesh methods into MeshOps.

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -27,36 +27,43 @@
 
 #include "pylith/utils/error.h" // USES PYLITH_CHECK_ERROR
 
-extern PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented);
+extern
+PetscErrorCode DMPlexGetOrientedFace(PetscDM dm,
+				     PetscInt cell,
+				     PetscInt faceSize,
+				     const PetscInt face[],
+				     PetscInt numCorners,
+				     PetscInt indices[],
+				     PetscInt origVertices[],
+				     PetscInt faceVertices[],
+				     PetscBool *posOriented);
 
 // ----------------------------------------------------------------------
 void
 pylith::faults::CohesiveTopology::createFault(topology::SubMesh* faultMesh,
-                                              DM& faultBoundary,
+                                              PetscDM& faultBoundary,
                                               const topology::Mesh& mesh,
-                                              DMLabel groupField,
+                                              PetscDMLabel groupField,
                                               const bool flipFault)
 { // createFault
   PYLITH_METHOD_BEGIN;
 
-  assert(0 != faultMesh);
+  assert(faultMesh);
   assert(groupField);
   PetscErrorCode err;
 
   faultMesh->coordsys(mesh);
-  DM       dmMesh = mesh.dmMesh();
-  PetscInt dim, depth;
+  PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
 
-  assert(dmMesh);
+  PetscInt dim, depth;
   err = DMPlexGetDimension(dmMesh, &dim);PYLITH_CHECK_ERROR(err);
   err = DMPlexGetDepth(dmMesh, &depth);PYLITH_CHECK_ERROR(err);
 
   // Convert fault to a DM
   if (depth == dim) {
-    DM                 subdm;
-    DMLabel            label;
-    const char        *groupName, *labelName = "boundary";
-    std::ostringstream tmp;
+    PetscDM subdm = NULL;
+    PetscDMLabel label = NULL;
+    const char *groupName = NULL, *labelName = "boundary";
 
     err = DMLabelGetName(groupField, &groupName);PYLITH_CHECK_ERROR(err);
     err = DMPlexCreateSubmesh(dmMesh, groupName, 1, &subdm);PYLITH_CHECK_ERROR(err);
@@ -64,14 +71,15 @@
     err = DMPlexGetLabel(subdm, labelName, &label);PYLITH_CHECK_ERROR(err);
     err = DMPlexMarkBoundaryFaces(subdm, label);PYLITH_CHECK_ERROR(err);
     err = DMPlexCreateSubmesh(subdm, labelName, 1, &faultBoundary);PYLITH_CHECK_ERROR(err);
-    faultMesh->setDMMesh(subdm);
+    std::string submeshLabel = "fault_" + std::string(groupName);
+    faultMesh->dmMesh(subdm, submeshLabel.c_str());
   } else {
-    DM              faultDMMeshTmp, faultDMMesh;
-    DMLabel         subpointMapTmp, subpointMap;
-    IS              pointIS;
-    const PetscInt *points;
-    PetscInt        depth, newDepth, h, numPoints, p;
-    const char     *groupName;
+    PetscDM faultDMMeshTmp = NULL, faultDMMesh = NULL;
+    PetscDMLabel subpointMapTmp = NULL, subpointMap = NULL;
+    PetscIS pointIS = NULL;
+    const PetscInt *points = NULL;
+    PetscInt depth, newDepth, h, numPoints, p;
+    const char *groupName = NULL;
 
     err = DMLabelGetName(groupField, &groupName);PYLITH_CHECK_ERROR(err);
     err = DMPlexCreateSubmesh(dmMesh, groupName, 1, &faultDMMeshTmp);PYLITH_CHECK_ERROR(err);
@@ -129,11 +137,13 @@
       err = DMRestoreWorkArray(faultDMMesh, maxConeSize, PETSC_INT, &revcone);PYLITH_CHECK_ERROR(err);
       err = DMRestoreWorkArray(faultDMMesh, maxConeSize, PETSC_INT, &revconeO);PYLITH_CHECK_ERROR(err);
     }
-    faultMesh->setDMMesh(faultDMMesh);
 
-    DMLabel            label;
-    const char        *labelName = "boundary";
+    std::string submeshLabel = "fault_" + std::string(groupName);
+    faultMesh->dmMesh(faultDMMesh, submeshLabel.c_str());
 
+    PetscDMLabel label = NULL;
+    const char *labelName = "boundary";
+
     err = DMPlexCreateLabel(faultDMMesh, labelName);PYLITH_CHECK_ERROR(err);
     err = DMPlexGetLabel(faultDMMesh, labelName, &label);PYLITH_CHECK_ERROR(err);
     err = DMPlexMarkBoundaryFaces(faultDMMesh, label);PYLITH_CHECK_ERROR(err);
@@ -157,7 +167,7 @@
 { // create
   PYLITH_METHOD_BEGIN;
 
-  assert(0 != mesh);
+  assert(mesh);
   assert(groupField);
 
   const char    *groupName;
@@ -614,90 +624,47 @@
   err = ISRestoreIndices(fVertexIS, &fVerticesDM);PYLITH_CHECK_ERROR(err);
   err = ISDestroy(&fVertexIS);PYLITH_CHECK_ERROR(err);
 
-  mesh->setDMMesh(newMesh);
+  mesh->dmMesh(newMesh);
 
   PYLITH_METHOD_END;
 } // create
 
 void
 pylith::faults::CohesiveTopology::createInterpolated(topology::Mesh* mesh,
-                                         const topology::SubMesh& faultMesh,
-                                         DM faultBoundary,
-                                         DMLabel groupField,
-                                         const int materialId,
-                                         int& firstFaultVertex,
-                                         int& firstLagrangeVertex,
-                                         int& firstFaultCell,
-                                         const bool constraintCell)
+						     const topology::SubMesh& faultMesh,
+						     PetscDM faultBoundary,
+						     PetscDMLabel groupField,
+						     const int materialId,
+						     int& firstFaultVertex,
+						     int& firstLagrangeVertex,
+						     int& firstFaultCell,
+						     const bool constraintCell)
 { // createInterpolated
-  assert(0 != mesh);
+  assert(mesh);
   assert(faultBoundary);
   assert(groupField);
-  DM             dm = mesh->dmMesh(), sdm;
-  DMLabel        label;
-  const char    *labelName = "faultSurface";
+  PetscDM sdm = NULL;
+  PetscDMLabel label = NULL;
+  const char *labelName = "faultSurface";
   PetscErrorCode err;
 
+  PetscDM dm = mesh->dmMesh();assert(dm);
   err = DMPlexGetLabel(dm, labelName, &label);PYLITH_CHECK_ERROR(err);
   // Completes the set of cells scheduled to be replaced
   //   Have to do internal fault vertices before fault boundary vertices, and this is the only thing I use faultBoundary for
   err = DMPlexLabelCohesiveComplete(dm, label);PYLITH_CHECK_ERROR(err);
   err = DMPlexConstructCohesiveCells(dm, label, &sdm);PYLITH_CHECK_ERROR(err);
-  mesh->setDMMesh(sdm);
+  mesh->dmMesh(sdm);
 } // createInterpolated
 
-PetscInt convertSieveToDMPointNumbering(PetscInt sievePoint, PetscInt numNormalCells, PetscInt numCohesiveCells, PetscInt numNormalVertices, PetscInt numShadowVertices, PetscInt numLagrangeVertices)
-{
-  PetscInt dmPoint = -1;
-
-  if ((sievePoint >= 0) && (sievePoint < numNormalCells)) {
-    dmPoint = sievePoint;
-    //std::cout << "normal cell sieve point "<<sievePoint<<" --> "<<" dm point"<<dmPoint<<std::endl;
-  } else if ((sievePoint >= numNormalCells) && (sievePoint < numNormalCells+numNormalVertices)) {
-    dmPoint = sievePoint+numCohesiveCells;
-    //std::cout << "normal vertex sieve point "<<sievePoint<<" --> "<<" dm point"<<dmPoint<<std::endl;
-  } else if ((sievePoint >= numNormalCells+numNormalVertices) && (sievePoint < numNormalCells+numNormalVertices+numShadowVertices+numLagrangeVertices)) {
-    dmPoint = sievePoint+numCohesiveCells;
-    //std::cout << "extra vertex sieve point "<<sievePoint<<" --> "<<" dm point"<<dmPoint<<std::endl;
-  } else if ((sievePoint >= numNormalCells+numNormalVertices+numShadowVertices+numLagrangeVertices) && (sievePoint < numNormalCells+numNormalVertices+numShadowVertices+numLagrangeVertices+numCohesiveCells)) {
-    dmPoint = sievePoint-(numNormalVertices+numShadowVertices+numLagrangeVertices);
-    //std::cout << "extra cell sieve point "<<sievePoint<<" --> "<<" dm point"<<dmPoint<<std::endl;
-  } else {
-    //std::cout << "face sieve point "<<sievePoint<<" --> "<<" dm point"<<dmPoint<<std::endl;
-  }
-  return dmPoint;
-}
-
-PetscInt convertDMToSievePointNumbering(PetscInt dmPoint, PetscInt numNormalCells, PetscInt numCohesiveCells, PetscInt numNormalVertices, PetscInt numShadowVertices, PetscInt numLagrangeVertices)
-{
-  PetscInt sievePoint = -1;
-
-  if ((dmPoint >= 0) && (dmPoint < numNormalCells)) {
-    sievePoint = dmPoint;
-    //std::cout << "normal cell sieve point "<<sievePoint<<" <-- "<<" dm point"<<dmPoint<<std::endl;
-  } else if ((dmPoint >= numNormalCells) && (dmPoint < numNormalCells+numCohesiveCells)) {
-    sievePoint = dmPoint+numNormalVertices+numShadowVertices+numLagrangeVertices;
-    //std::cout << "extra cell sieve point "<<sievePoint<<" <-- "<<" dm point"<<dmPoint<<std::endl;
-  } else if ((dmPoint >= numNormalCells+numCohesiveCells) && (dmPoint < numNormalCells+numCohesiveCells+numNormalVertices)) {
-    sievePoint = dmPoint-numCohesiveCells;
-    //std::cout << "normal vertex sieve point "<<sievePoint<<" <-- "<<" dm point"<<dmPoint<<std::endl;
-  } else if ((dmPoint >= numNormalCells+numCohesiveCells+numNormalVertices) && (dmPoint < numNormalCells+numCohesiveCells+numNormalVertices+numShadowVertices+numLagrangeVertices)) {
-    sievePoint = dmPoint-numCohesiveCells;
-    //std::cout << "extra vertex sieve point "<<sievePoint<<" <-- "<<" dm point"<<dmPoint<<std::endl;
-  } else {
-    //std::cout << "face sieve point "<<sievePoint<<" <-- "<<" dm point"<<dmPoint<<std::endl;
-  }
-  return sievePoint;
-}
-
 // ----------------------------------------------------------------------
 // Form a parallel fault mesh using the cohesive cell information
 void
-pylith::faults::CohesiveTopology::createFaultParallel(
-			    topology::SubMesh* faultMesh,
-			    const topology::Mesh& mesh,
-			    const int materialId,
-			    const bool constraintCell)
+pylith::faults::CohesiveTopology::createFaultParallel(topology::SubMesh* faultMesh,
+						      const topology::Mesh& mesh,
+						      const int materialId,
+						      const char* label,
+						      const bool constraintCell)
 { // createFaultParallel
   PYLITH_METHOD_BEGIN;
 
@@ -706,12 +673,12 @@
 
   faultMesh->coordsys(mesh);
 
-  DM dmMesh = mesh.dmMesh();
-  assert(dmMesh);
-  DM dmFaultMesh;
+  PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
+  PetscDM dmFaultMesh;
 
   err = DMPlexCreateCohesiveSubmesh(dmMesh, constraintCell ? PETSC_TRUE : PETSC_FALSE, &dmFaultMesh);PYLITH_CHECK_ERROR(err);
-  faultMesh->setDMMesh(dmFaultMesh);
+  std::string meshLabel = "fault_" + std::string(label);
+  faultMesh->dmMesh(dmFaultMesh, meshLabel.c_str());
 
   PYLITH_METHOD_END;
 } // createFaultParallel

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.hh	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.hh	2013-05-13 18:43:34 UTC (rev 22057)
@@ -112,6 +112,7 @@
    * @param faultMesh Finite-element mesh of fault (output).
    * @param mesh Finite-element mesh.
    * @param materialId Material id for cohesive elements.
+   * @param label Fault label.
    * @param constraintCell True if creating cells constrained with 
    *   Lagrange multipliers that require extra vertices, false otherwise.
    */
@@ -119,6 +120,7 @@
   void createFaultParallel(topology::SubMesh* faultMesh,
 			   const topology::Mesh& mesh,
 			   const int materialId,
+			   const char* label,
 			   const bool constraintCell =false);
 
 }; // class CohesiveTopology

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -92,7 +92,7 @@
   const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();assert(cs);
 
   delete _faultMesh; _faultMesh = new topology::SubMesh();assert(_faultMesh);
-  CohesiveTopology::createFaultParallel(_faultMesh, mesh, id(), _useLagrangeConstraints);
+  CohesiveTopology::createFaultParallel(_faultMesh, mesh, id(), label(), _useLagrangeConstraints);
   _initializeCohesiveInfo(mesh);
 
   delete _fields;

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -67,8 +67,7 @@
   assert(_quadrature);
 
   delete _faultMesh; _faultMesh = new topology::SubMesh();
-  CohesiveTopology::createFaultParallel(_faultMesh, mesh, id(),
-    useLagrangeConstraints());
+  CohesiveTopology::createFaultParallel(_faultMesh, mesh, id(), label(), useLagrangeConstraints());
 
   // Reset fields.
   delete _fields; 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -77,7 +77,7 @@
 
   err = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);PYLITH_CHECK_ERROR(err);
   err = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, pInterpolate, &cells[0], spaceDim, &(*coordinates)[0], &dmMesh);PYLITH_CHECK_ERROR(err);
-  mesh->setDMMesh(dmMesh);
+  mesh->dmMesh(dmMesh);
 
   PYLITH_METHOD_END;
 } // buildMesh

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -21,6 +21,7 @@
 #include "OutputSolnPoints.hh" // implementation of class methods
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
 #include "MeshBuilder.hh" // USES MeshBuilder
@@ -125,7 +126,7 @@
   // Create mesh corresponding to points.
   const int meshDim = 0;
   delete _pointsMesh; _pointsMesh = new topology::Mesh(meshDim);assert(_pointsMesh);
-  _pointsMesh->createDMMesh(meshDim);
+  topology::MeshOps::createDMMesh(_pointsMesh, meshDim, _mesh->comm(), "points");
 
   const int numPointsLocal = _interpolator->n;
   PylithScalar* pointsLocal = NULL;
@@ -150,7 +151,7 @@
 
   // Set coordinate system and create nondimensionalized coordinates
   _pointsMesh->coordsys(_mesh->coordsys());
-  _pointsMesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(_pointsMesh, normalizer);
   
 #if 0 // DEBUGGING
   _pointsMesh->view("POINTS MESH");

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -61,7 +61,7 @@
 
   PetscDM newDM = NULL;
   PetscErrorCode err = DMPlexDistribute(origMesh.dmMesh(), partitioner, 0, &newDM);PYLITH_CHECK_ERROR(err);
-  newMesh->setDMMesh(newDM);
+  newMesh->dmMesh(newDM);
 
   PYLITH_METHOD_END;
 } // distribute

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -22,7 +22,6 @@
 #include "Mesh.hh" // implementation of class methods
 
 #include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
-#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 #include "pylith/utils/array.hh" // USES scalar_array
 #include "pylith/utils/petscfwd.h" // USES PetscVec
 #include "pylith/utils/error.h" // USES PYLITH_CHECK_ERROR
@@ -34,31 +33,106 @@
 // ----------------------------------------------------------------------
 // Default constructor
 pylith::topology::Mesh::Mesh(void) :
-  _newMesh(NULL),
-  _numNormalCells(0), _numCohesiveCells(0), _numNormalVertices(0), _numShadowVertices(0), _numLagrangeVertices(0),
+  _dmMesh(NULL),
+  _numNormalCells(0),
+  _numCohesiveCells(0),
+  _numNormalVertices(0),
+  _numShadowVertices(0),
+  _numLagrangeVertices(0),
   _coordsys(0),
   _comm(PETSC_COMM_WORLD),
-  _debug(false)
+  _debug(false),
+  _isSubMesh(false)
 { // constructor
 } // constructor
 
 // ----------------------------------------------------------------------
-// Default constructor
+// Constructor with dimension and communicator.
 pylith::topology::Mesh::Mesh(const int dim,
 			     const MPI_Comm& comm) :
-  _newMesh(NULL),
-  _numNormalCells(0), _numCohesiveCells(0), _numNormalVertices(0), _numShadowVertices(0), _numLagrangeVertices(0),
+  _dmMesh(NULL),
+  _numNormalCells(0),
+  _numCohesiveCells(0),
+  _numNormalVertices(0),
+  _numShadowVertices(0),
+  _numLagrangeVertices(0),
   _coordsys(0),
   _comm(comm),
-  _debug(false)
+  _debug(false),
+  _isSubMesh(false)
 { // constructor
   PYLITH_METHOD_BEGIN;
 
-  createDMMesh(dim);
+  PetscErrorCode err;
+  err = DMCreate(_comm, &_dmMesh);PYLITH_CHECK_ERROR(err);
+  err = DMSetType(_dmMesh, DMPLEX);PYLITH_CHECK_ERROR(err);
+  err = DMPlexSetDimension(_dmMesh, dim);PYLITH_CHECK_ERROR(err);
+  err = PetscObjectSetName((PetscObject) _dmMesh, "domain");PYLITH_CHECK_ERROR(err);
+
   PYLITH_METHOD_END;
 } // constructor
 
 // ----------------------------------------------------------------------
+// Create submesh.
+pylith::topology::Mesh::Mesh(const Mesh& mesh,
+			     const char* label) :
+  _dmMesh(NULL),
+  _numNormalCells(0),
+  _numCohesiveCells(0),
+  _numNormalVertices(0),
+  _numShadowVertices(0),
+  _numLagrangeVertices(0),
+  _coordsys(0),
+  _comm(mesh.comm()),
+  _debug(mesh.debug()),
+  _isSubMesh(true)
+{ // Submesh constructor
+  PYLITH_METHOD_BEGIN;
+
+  assert(label);
+
+  coordsys(mesh.coordsys());
+
+  PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
+  PetscErrorCode err;
+
+  PetscBool hasLabel = PETSC_FALSE;
+  err = DMPlexHasLabel(dmMesh, label, &hasLabel);PYLITH_CHECK_ERROR(err);
+  if (!hasLabel) {
+    std::ostringstream msg;
+    msg << "Could not find group of points '" << label << "' in PETSc DM mesh.";
+    throw std::runtime_error(msg.str());
+  } // if
+
+  if (mesh.dimension() < 1) {
+    throw std::logic_error("INTERNAL ERROR in MeshOps::createSubMesh()\n"
+			   "Cannot create submesh for mesh with dimension < 1.");
+  } // if
+
+  /* TODO: Add creation of pointSF for submesh */
+  err = DMPlexCreateSubmesh(dmMesh, label, 1, &_dmMesh);PYLITH_CHECK_ERROR(err);
+
+  PetscInt maxConeSizeLocal = 0, maxConeSize = 0;
+  err = DMPlexGetMaxSizes(_dmMesh, &maxConeSizeLocal, NULL);PYLITH_CHECK_ERROR(err);
+  err = MPI_Allreduce(&maxConeSizeLocal, &maxConeSize, 1, MPI_INT, MPI_MAX,
+                      PetscObjectComm((PetscObject) _dmMesh)); PYLITH_CHECK_ERROR(err);
+
+  if (maxConeSize <= 0) {
+    std::ostringstream msg;
+    msg << "Error while creating submesh. Submesh '" 
+	<< label << "' does not contain any cells.\n";
+    throw std::runtime_error(msg.str());
+  } // if
+
+  // Set name
+  std::string meshLabel = "subdomain_" + std::string(label);
+  err = PetscObjectSetName((PetscObject) _dmMesh, meshLabel.c_str());PYLITH_CHECK_ERROR(err);
+
+  PYLITH_METHOD_END;
+} // SubMesh constructor
+		     
+
+// ----------------------------------------------------------------------
 // Default destructor
 pylith::topology::Mesh::~Mesh(void)
 { // destructor
@@ -73,29 +147,12 @@
   PYLITH_METHOD_BEGIN;
 
   delete _coordsys; _coordsys = 0;
-  PetscErrorCode err = DMDestroy(&_newMesh);PYLITH_CHECK_ERROR(err);
+  PetscErrorCode err = DMDestroy(&_dmMesh);PYLITH_CHECK_ERROR(err);
 
   PYLITH_METHOD_END;
 } // deallocate
   
 // ----------------------------------------------------------------------
-// Create DMPlex mesh.
-void
-pylith::topology::Mesh::createDMMesh(const int dim)
-{ // createDMMesh
-  PYLITH_METHOD_BEGIN;
-
-  PetscErrorCode err;
-  err = DMDestroy(&_newMesh);PYLITH_CHECK_ERROR(err);
-  err = DMCreate(_comm, &_newMesh);PYLITH_CHECK_ERROR(err);
-  err = DMSetType(_newMesh, DMPLEX);PYLITH_CHECK_ERROR(err);
-  err = DMPlexSetDimension(_newMesh, dim);PYLITH_CHECK_ERROR(err);
-  err = PetscObjectSetName((PetscObject) _newMesh, "domain");PYLITH_CHECK_ERROR(err);
-
-  PYLITH_METHOD_END;
-} // createDMMesh
-
-// ----------------------------------------------------------------------
 // Set coordinate system.
 void
 pylith::topology::Mesh::coordsys(const spatialdata::geocoords::CoordSys* cs)
@@ -122,17 +179,17 @@
   *numNames = 0;
   *names = 0;
 
-  if (_newMesh) {
+  if (_dmMesh) {
     PetscErrorCode err = 0;
 
     PetscInt numLabels = 0;
-    err = DMPlexGetNumLabels(_newMesh, &numLabels);PYLITH_CHECK_ERROR(err);
+    err = DMPlexGetNumLabels(_dmMesh, &numLabels);PYLITH_CHECK_ERROR(err);
 
     *numNames = numLabels;
     *names = new char*[numLabels];
     for (int iLabel=0; iLabel < numLabels; ++iLabel) {
       const char* namestr = NULL;
-      err = DMPlexGetLabelName(_newMesh, iLabel, &namestr);PYLITH_CHECK_ERROR(err);
+      err = DMPlexGetLabelName(_dmMesh, iLabel, &namestr);PYLITH_CHECK_ERROR(err);
       // Must return char* that SWIG can deallocate.
       const char len = strlen(namestr);
       char* newName = 0;
@@ -157,12 +214,12 @@
 { // groupSize
   PYLITH_METHOD_BEGIN;
 
-  assert(_newMesh);
+  assert(_dmMesh);
 
   PetscErrorCode err = 0;
 
   PetscBool hasLabel = PETSC_FALSE;
-  err = DMPlexHasLabel(_newMesh, name, &hasLabel);PYLITH_CHECK_ERROR(err);
+  err = DMPlexHasLabel(_dmMesh, name, &hasLabel);PYLITH_CHECK_ERROR(err);
   if (!hasLabel) {
     std::ostringstream msg;
     msg << "Cannot get size of group '" << name
@@ -171,31 +228,10 @@
   } // if
 
   PetscInt size = 0;
-  err = DMPlexGetLabelSize(_newMesh, name, &size);PYLITH_CHECK_ERROR(err);
+  err = DMPlexGetLabelSize(_dmMesh, name, &size);PYLITH_CHECK_ERROR(err);
 
   PYLITH_METHOD_RETURN(size);
 } // groupSize
 
 
-// ----------------------------------------------------------------------
-// Nondimensionalize the finite-element mesh.
-void 
-pylith::topology::Mesh::nondimensionalize(const spatialdata::units::Nondimensional& normalizer)
-{ // initialize
-  PYLITH_METHOD_BEGIN;
-
-  PetscVec coordVec, coordDimVec;
-  const PylithScalar lengthScale = normalizer.lengthScale();
-  PetscErrorCode err;
-
-  assert(_newMesh);
-  err = DMGetCoordinatesLocal(_newMesh, &coordVec);PYLITH_CHECK_ERROR(err);assert(coordVec);
-  // There does not seem to be an advantage to calling nondimensionalize()
-  err = VecScale(coordVec, 1.0/lengthScale);PYLITH_CHECK_ERROR(err);
-  err = DMPlexSetScale(_newMesh, PETSC_UNIT_LENGTH, lengthScale);PYLITH_CHECK_ERROR(err);
-
-  PYLITH_METHOD_END;
-} // nondimensionalize
-
-
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.hh	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.hh	2013-05-13 18:43:34 UTC (rev 22057)
@@ -28,7 +28,6 @@
 // Include directives ---------------------------------------------------
 #include "topologyfwd.hh" // forward declarations
 #include "spatialdata/geocoords/geocoordsfwd.hh" // forward declarations
-#include "spatialdata/units/unitsfwd.hh" // forward declarations
 
 #include "pylith/utils/petscfwd.h" // HASA PetscDM
 
@@ -56,18 +55,20 @@
   Mesh(const int dim,
        const MPI_Comm& comm =PETSC_COMM_WORLD); 
 
+  /** Create submesh.
+   *
+   * @param mesh Mesh over domain.
+   * @param label Label of vertices on boundary.
+   */
+  Mesh(const Mesh& mesh,
+       const char* label);
+
   /// Default destructor
   ~Mesh(void);
 
   /// Deallocate PETSc and local data structures.
   void deallocate(void);
   
-  /** Create DMPlex mesh.
-   *
-   * @param dim Dimension associated with mesh cells.
-   */
-  void createDMMesh(const int dim=3); 
-
   /** Get DMPlex mesh.
    *
    * @returns DMPlex mesh.
@@ -77,8 +78,10 @@
   /** Set DMPlex mesh.
    *
    * @param DMPlex mesh.
+   * @param label Label for mesh.
    */
-  void setDMMesh(PetscDM dm);
+  void dmMesh(PetscDM dm,
+	      const char* label ="domain");
 
   /** Get sizes for all point types.
    *
@@ -88,7 +91,11 @@
    * @param numShadowVertices
    * @param numLagrangeVertices.
    */
-  void getPointTypeSizes(PetscInt *numNormalCells, PetscInt *numCohesiveCells, PetscInt *numNormalVertices, PetscInt *numShadowVertices, PetscInt *numLagrangeVertices) const;
+  void getPointTypeSizes(PetscInt *numNormalCells, 
+			 PetscInt *numCohesiveCells,
+			 PetscInt *numNormalVertices,
+			 PetscInt *numShadowVertices,
+			 PetscInt *numLagrangeVertices) const;
 
   /** Set sizes for all point types.
    *
@@ -98,7 +105,11 @@
    * @param numShadowVertices
    * @param numLagrangeVertices.
    */
-  void setPointTypeSizes(PetscInt numNormalCells, PetscInt numCohesiveCells, PetscInt numNormalVertices, PetscInt numShadowVertices, PetscInt numLagrangeVertices);
+  void setPointTypeSizes(PetscInt numNormalCells,
+			 PetscInt numCohesiveCells,
+			 PetscInt numNormalVertices,
+			 PetscInt numShadowVertices,
+			 PetscInt numLagrangeVertices);
 
   /** Set coordinate system.
    *
@@ -130,18 +141,18 @@
    */
   int dimension(void) const;
 
+  /** Get number of vertices in mesh.
+   *
+   * @returns Number of vertices in mesh.
+   */
+  int numVertices(void) const;
+  
   /** Get representative cone size for mesh.
    *
    * @returns Representative cone size for mesh.
    */
   int coneSize(void) const;
   
-  /** Get number of vertices in mesh.
-   *
-   * @returns Number of vertices in mesh.
-   */
-  int numVertices(void) const;
-  
   /** Get number of cells in mesh.
    *
    * @returns Number of cells in mesh.
@@ -188,16 +199,10 @@
    */
   int groupSize(const char *name);
 
-  /** Nondimensionalize the finite-element mesh.
-   *
-   * @param normalizer Nondimensionalizer.
-   */
-  void nondimensionalize(const spatialdata::units::Nondimensional& normalizer);
-
 // PRIVATE MEMBERS //////////////////////////////////////////////////////
 private :
 
-  PetscDM _newMesh;
+  PetscDM _dmMesh;
 
   /* The old-style point numbering: [normalCells, normalVertices, shadowVertices, lagrangeVertices, cohesiveCells]
      The new-style point numbering: [normalCells, cohesiveCells, normalVertices, shadowVertices, lagrangeVertices]
@@ -207,6 +212,7 @@
   spatialdata::geocoords::CoordSys* _coordsys; ///< Coordinate system.
   MPI_Comm _comm; ///< MPI communicator for mesh.
   bool _debug; ///< Debugging flag for mesh.
+  bool _isSubMesh; ///< True if mesh is a submesh of another mesh.
   
 // NOT IMPLEMENTED //////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -22,49 +22,63 @@
 
 #include "pylith/utils/error.h" // USES PYLITH_CHECK_ERROR
 
+// ----------------------------------------------------------------------
 // Get DMPlex mesh.
 inline
-DM
+PetscDM
 pylith::topology::Mesh::dmMesh(void) const {
-  return _newMesh;
+  return _dmMesh;
 }
 
+// ----------------------------------------------------------------------
 // Set DMPlex mesh.
 inline
 void
-pylith::topology::Mesh::setDMMesh(DM dm) {
+pylith::topology::Mesh::dmMesh(PetscDM dm,
+			       const char* label) {
   PYLITH_METHOD_BEGIN;
 
   PetscErrorCode err;
-  err = DMDestroy(&_newMesh);PYLITH_CHECK_ERROR(err);
-  _newMesh = dm;
-  err = PetscObjectSetName((PetscObject) _newMesh, "domain");PYLITH_CHECK_ERROR(err);
+  err = DMDestroy(&_dmMesh);PYLITH_CHECK_ERROR(err);
+  _dmMesh = dm;
+  err = PetscObjectSetName((PetscObject) _dmMesh, label);PYLITH_CHECK_ERROR(err);
 
   PYLITH_METHOD_END;
 }
 
+// ----------------------------------------------------------------------
 // Get point type sizes.
 inline
 void
-pylith::topology::Mesh::getPointTypeSizes(PetscInt *numNormalCells, PetscInt *numCohesiveCells, PetscInt *numNormalVertices, PetscInt *numShadowVertices, PetscInt *numLagrangeVertices) const {
-  *numNormalCells      = _numNormalCells;
-  *numCohesiveCells    = _numCohesiveCells;
-  *numNormalVertices   = _numNormalVertices;
-  *numShadowVertices   = _numShadowVertices;
+pylith::topology::Mesh::getPointTypeSizes(PetscInt *numNormalCells,
+					  PetscInt *numCohesiveCells,
+					  PetscInt *numNormalVertices,
+					  PetscInt *numShadowVertices,
+					  PetscInt *numLagrangeVertices) const {
+  *numNormalCells = _numNormalCells;
+  *numCohesiveCells = _numCohesiveCells;
+  *numNormalVertices = _numNormalVertices;
+  *numShadowVertices = _numShadowVertices;
   *numLagrangeVertices = _numLagrangeVertices;
 }
 
+// ----------------------------------------------------------------------
 // Set point type sizes.
 inline
 void
-pylith::topology::Mesh::setPointTypeSizes(PetscInt numNormalCells, PetscInt numCohesiveCells, PetscInt numNormalVertices, PetscInt numShadowVertices, PetscInt numLagrangeVertices) {
-  _numNormalCells      = numNormalCells;
-  _numCohesiveCells    = numCohesiveCells;
-  _numNormalVertices   = numNormalVertices;
-  _numShadowVertices   = numShadowVertices;
+pylith::topology::Mesh::setPointTypeSizes(PetscInt numNormalCells,
+					  PetscInt numCohesiveCells,
+					  PetscInt numNormalVertices,
+					  PetscInt numShadowVertices,
+					  PetscInt numLagrangeVertices) {
+  _numNormalCells = numNormalCells;
+  _numCohesiveCells = numCohesiveCells;
+  _numNormalVertices = numNormalVertices;
+  _numShadowVertices = numShadowVertices;
   _numLagrangeVertices = numLagrangeVertices;
 }
 
+// ----------------------------------------------------------------------
 // Get coordinate system.
 inline
 const spatialdata::geocoords::CoordSys*
@@ -72,6 +86,7 @@
   return _coordsys;
 }
 
+// ----------------------------------------------------------------------
 // Set debug flag.
 inline
 void
@@ -79,6 +94,7 @@
   _debug = value;
 }
 
+// ----------------------------------------------------------------------
 // Get debug flag.
 inline
 bool
@@ -86,6 +102,7 @@
   return _debug;
 }
 
+// ----------------------------------------------------------------------
 // Get dimension of mesh.
 inline
 int
@@ -93,46 +110,50 @@
   PYLITH_METHOD_BEGIN;
 
   PetscInt dim = 0;
-  if (_newMesh) {
-    PetscErrorCode err = DMPlexGetDimension(_newMesh, &dim);PYLITH_CHECK_ERROR(err);
+  if (_dmMesh) {
+    PetscErrorCode err = DMPlexGetDimension(_dmMesh, &dim);PYLITH_CHECK_ERROR(err);
   } // if
 
   PYLITH_METHOD_RETURN(dim);
 }
 
-// Get representative cone size for mesh.
+// ----------------------------------------------------------------------
+// Get number of vertices in mesh.
 inline
 int
-pylith::topology::Mesh::coneSize(void) const {
+pylith::topology::Mesh::numVertices(void) const {
   PYLITH_METHOD_BEGIN;
 
-  PetscInt coneSize = 0;
-  if (_newMesh) {
-    PetscInt cStart;
-    PetscErrorCode err;
-    err = DMPlexGetHeightStratum(_newMesh, 0, &cStart, NULL);PYLITH_CHECK_ERROR(err);
-    err = DMPlexGetConeSize(_newMesh, cStart, &coneSize);PYLITH_CHECK_ERROR(err);
+  PetscInt nvertices = 0;
+  if (_dmMesh) {
+    PetscInt begin=0, end=0;
+    PetscErrorCode err = DMPlexGetDepthStratum(_dmMesh, 0, &begin, &end);PYLITH_CHECK_ERROR(err);
+    nvertices = end-begin;
   } // if
 
-  PYLITH_METHOD_RETURN(coneSize);
+  PYLITH_METHOD_RETURN(nvertices);
 }
 
-// Get number of vertices in mesh.
+// ----------------------------------------------------------------------
+// Get representative cone size for mesh.
 inline
 int
-pylith::topology::Mesh::numVertices(void) const {
+pylith::topology::Mesh::coneSize(void) const {
   PYLITH_METHOD_BEGIN;
 
-  PetscInt nvertices = 0;
-  if (_newMesh) {
-    PetscInt begin=0, end=0;
-    PetscErrorCode err = DMPlexGetDepthStratum(_newMesh, 0, &begin, &end);PYLITH_CHECK_ERROR(err);
-    nvertices = end-begin;
+  PetscInt coneSize = 0;
+  if (_dmMesh) {
+    PetscInt cStart;
+    PetscErrorCode err;
+    const int cellHeight = _isSubMesh ? 1 : 0;
+    err = DMPlexGetHeightStratum(_dmMesh, cellHeight, &cStart, NULL);PYLITH_CHECK_ERROR(err);
+    err = DMPlexGetConeSize(_dmMesh, cStart, &coneSize);PYLITH_CHECK_ERROR(err);
   } // if
 
-  PYLITH_METHOD_RETURN(nvertices);
+  PYLITH_METHOD_RETURN(coneSize);
 }
 
+// ----------------------------------------------------------------------
 // Get number of cells in mesh.
 inline
 int
@@ -140,15 +161,17 @@
   PYLITH_METHOD_BEGIN;
 
   PetscInt ncells = 0;
-  if (_newMesh) {
+  if (_dmMesh) {
     PetscInt begin=0, end=0;
-    PetscErrorCode err = DMPlexGetHeightStratum(_newMesh, 0, &begin, &end);PYLITH_CHECK_ERROR(err);
+    const int cellHeight = _isSubMesh ? 1 : 0;
+    PetscErrorCode err = DMPlexGetHeightStratum(_dmMesh, cellHeight, &begin, &end);PYLITH_CHECK_ERROR(err);
     ncells = end-begin;
   } // if
 
   PYLITH_METHOD_RETURN(ncells);
 }
 
+// ----------------------------------------------------------------------
 // Set MPI communicator associated with mesh.
 inline
 void
@@ -156,6 +179,7 @@
   _comm = value;
 }
     
+// ----------------------------------------------------------------------
 // Get MPI communicator associated with mesh.
 inline
 const MPI_Comm
@@ -163,6 +187,7 @@
   return _comm;
 }
     
+// ----------------------------------------------------------------------
 // Get MPI rank.
 inline
 int
@@ -172,14 +197,15 @@
   return rank;
 }
     
+// ----------------------------------------------------------------------
 // Print mesh to stdout.
 inline
 void
 pylith::topology::Mesh::view(const char* label) const {
   PYLITH_METHOD_BEGIN;
 
-  assert(_newMesh);
-  PetscErrorCode err = DMView(_newMesh, PETSC_VIEWER_STDOUT_WORLD);PYLITH_CHECK_ERROR(err);
+  assert(_dmMesh);
+  PetscErrorCode err = DMView(_dmMesh, PETSC_VIEWER_STDOUT_WORLD);PYLITH_CHECK_ERROR(err);
 
   PYLITH_METHOD_END;
 }

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -24,6 +24,8 @@
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/utils/array.hh" // USES int_array
 
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
 #include <stdexcept> // USES std::runtime_error
 #include <sstream> // USES std::ostringstream
 #include <cassert> // USES assert()
@@ -33,7 +35,50 @@
 
 
 // ----------------------------------------------------------------------
+// Create PETSc DM mesh.
 void
+pylith::topology::MeshOps::createDMMesh(Mesh* const mesh,
+					const int dim,
+					const MPI_Comm& comm,
+					const char* label)
+{ // createDMMesh
+  PYLITH_METHOD_BEGIN;
+
+  PetscErrorCode err;
+  PetscDM dmMesh = NULL;
+  err = DMCreate(comm, &dmMesh);PYLITH_CHECK_ERROR(err);
+  err = DMSetType(dmMesh, DMPLEX);PYLITH_CHECK_ERROR(err);
+  err = DMPlexSetDimension(dmMesh, dim);PYLITH_CHECK_ERROR(err);
+  mesh->dmMesh(dmMesh, label);
+  mesh->comm(comm);
+
+  PYLITH_METHOD_END;
+} // createDMMesh
+
+// ----------------------------------------------------------------------
+// Nondimensionalize the finite-element mesh.
+void 
+pylith::topology::MeshOps::nondimensionalize(Mesh* const mesh,
+					     const spatialdata::units::Nondimensional& normalizer)
+{ // nondimensionalize
+  PYLITH_METHOD_BEGIN;
+
+  PetscVec coordVec, coordDimVec;
+  const PylithScalar lengthScale = normalizer.lengthScale();
+  PetscErrorCode err;
+
+  PetscDM dmMesh = mesh->dmMesh();assert(dmMesh);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);PYLITH_CHECK_ERROR(err);assert(coordVec);
+  // There does not seem to be an advantage to calling nondimensionalize()
+  err = VecScale(coordVec, 1.0/lengthScale);PYLITH_CHECK_ERROR(err);
+  err = DMPlexSetScale(dmMesh, PETSC_UNIT_LENGTH, lengthScale);PYLITH_CHECK_ERROR(err);
+
+  PYLITH_METHOD_END;
+} // nondimensionalize
+
+
+// ----------------------------------------------------------------------
+void
 pylith::topology::MeshOps::checkMaterialIds(const Mesh& mesh,
 					    int* const materialIds,
 					    const int numMaterials)

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.hh	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.hh	2013-05-13 18:43:34 UTC (rev 22057)
@@ -28,6 +28,8 @@
 // Include directives ---------------------------------------------------
 #include "topologyfwd.hh" // forward declarations
 
+#include "spatialdata/units/unitsfwd.hh" // forward declarations
+
 // MeshOps --------------------------------------------------------------
 /// Simple operations on a Mesh object.
 class pylith::topology::MeshOps
@@ -37,6 +39,28 @@
 // PUBLIC MEMBERS ///////////////////////////////////////////////////////
 public :
 
+  /** Create DMPlex mesh.
+   *
+   * @param mesh Finite-element mesh.
+   * @param dim Dimension associated with mesh cells.
+   * @param comm MPI communicator for mesh.
+   * @param label Label for mesh.
+   */
+  static
+  void createDMMesh(Mesh* const mesh,
+		    const int dim =3,
+		    const MPI_Comm& comm =PETSC_COMM_WORLD,
+		    const char* label ="domain");
+
+  /** Nondimensionalize the finite-element mesh.
+   *
+   * @param mesh Finite-element mesh.
+   * @param normalizer Nondimensionalizer.
+   */
+  static
+  void nondimensionalize(Mesh* const mesh,
+			 const spatialdata::units::Nondimensional& normalizer);
+
   /** Check to make sure material id of every cell matches the id of
    *  one of the materials.
    *

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/ReverseCuthillMcKee.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/ReverseCuthillMcKee.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/ReverseCuthillMcKee.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -51,7 +51,7 @@
     sieveMesh->relabel(*reordering);
     //sieveMesh->view("MESH AFTER RELABEL");
 #else
-    assert(0);
+    std::cerr << "WARNING: Reverse Cuthill McKee temporarily disabled." << std::endl;
 #endif
   } // if    
 } // reorder

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -144,6 +144,7 @@
 { // initialize
   PYLITH_METHOD_BEGIN;
 
+  assert(0);
   if (_coordsys)
     _coordsys->initialize();
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh	2013-05-13 18:43:34 UTC (rev 22057)
@@ -79,8 +79,10 @@
   /** Set DMPlex mesh.
    *
    * @param DMPlex mesh.
+   * @param label Label for mesh.
    */
-  void setDMMesh(PetscDM dm);
+  void dmMesh(PetscDM dm,
+	      const char* label ="domain");
 
   /** Get sizes for all point types.
    *

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -32,11 +32,12 @@
 // Set DMComplex mesh.
 inline
 void
-pylith::topology::SubMesh::setDMMesh(PetscDM dm) {
+pylith::topology::SubMesh::dmMesh(PetscDM dm,
+				  const char* label) {
   PetscErrorCode err;
   err = DMDestroy(&_newMesh);PYLITH_CHECK_ERROR(err);
   _newMesh = dm;
-  err = PetscObjectSetName((PetscObject) _newMesh, "domain");PYLITH_CHECK_ERROR(err);
+  err = PetscObjectSetName((PetscObject) _newMesh, label);PYLITH_CHECK_ERROR(err);
 }
 
 // Get point type sizes.

Modified: short/3D/PyLith/trunk/modulesrc/topology/Mesh.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/Mesh.i	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/modulesrc/topology/Mesh.i	2013-05-13 18:43:34 UTC (rev 22057)
@@ -43,6 +43,14 @@
       Mesh(const int dim,
 	   const MPI_Comm& comm =PETSC_COMM_WORLD); 
 
+      /** Create submesh.
+       *
+       * @param mesh Mesh over domain.
+       * @param label Label of vertices on boundary.
+       */
+      Mesh(const Mesh& mesh,
+	   const char* label);
+
       /// Default destructor
       ~Mesh(void);
       
@@ -50,12 +58,6 @@
       virtual
       void deallocate(void);
 
-      /** Create DM mesh.
-       *
-       * @param dim Dimension associated with mesh cells.
-       */
-      void createDMMesh(const int dim=3); 
-      
       /** Set coordinate system.
        *
        * @param cs Coordinate system.
@@ -116,12 +118,6 @@
        */
       const MPI_Comm comm(void) const;
     
-      /** Initialize the finite-element mesh.
-       *
-       * @param normalizer Nondimensionalizer.
-       */
-      void nondimensionalize(const spatialdata::units::Nondimensional& normalizer);
-
       /** Print mesh to stdout.
        *
        * @param label Label for mesh.

Modified: short/3D/PyLith/trunk/modulesrc/topology/MeshOps.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/MeshOps.i	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/modulesrc/topology/MeshOps.i	2013-05-13 18:43:34 UTC (rev 22057)
@@ -22,6 +22,19 @@
  * @brief Python interface to C++ MeshOps.
  */
 
+%inline %{
+  /** Nondimensionalize the finite-element mesh.
+   *
+   * @param mesh Finite-element mesh.
+   * @param normalizer Nondimensionalizer.
+   */
+  void
+  MeshOps_nondimensionalize(pylith::topology::Mesh* const mesh,
+			    const spatialdata::units::Nondimensional& normalizer) {
+    pylith::topology::MeshOps::nondimensionalize(mesh, normalizer);
+  } // nondimensionalize
+%}
+
 %apply(int* IN_ARRAY1, int DIM1) {
   (int* const materialIds, const int numMaterials)
   };

Modified: short/3D/PyLith/trunk/pylith/topology/MeshImporter.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/MeshImporter.py	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/pylith/topology/MeshImporter.py	2013-05-13 18:43:34 UTC (rev 22057)
@@ -135,7 +135,8 @@
       newMesh.memLoggingStage = "RefinedMesh"
 
     # Nondimensionalize mesh (coordinates of vertices).
-    newMesh.nondimensionalize(normalizer)
+    from pylith.topology.topology import MeshOps_nondimensionalize
+    MeshOps_nondimensionalize(newMesh, normalizer)
 
     self._eventLogger.eventEnd(logEvent)    
     return newMesh

Modified: short/3D/PyLith/trunk/pylith/topology/MeshImporterDist.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/MeshImporterDist.py	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/pylith/topology/MeshImporterDist.py	2013-05-13 18:43:34 UTC (rev 22057)
@@ -94,7 +94,8 @@
     mesh = self.refiner.refine(mesh)
 
     # Nondimensionalize mesh (coordinates of vertices).
-    mesh.nondimensionalize(normalizer)
+    from pylith.topology.topology import MeshOps_nondimensionalize
+    MeshOps_nondimensionalize(mesh, normalizer)
 
     self._eventLogger.eventEnd(logEvent)    
     return mesh

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -25,6 +25,7 @@
 #include "data/AbsorbingDampersData.hh" // USES AbsorbingDampersData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Fields.hh" // USES Fields
@@ -421,7 +422,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Set up quadrature
   _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryConditionPoints.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryConditionPoints.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryConditionPoints.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -25,6 +25,7 @@
 #include "data/PointForceDataTri3.hh" // USES PointForceDataTri3
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 
@@ -53,7 +54,7 @@
   cs.setSpaceDim(mesh.dimension());
   cs.initialize();
   mesh.coordsys(&cs);
-  mesh.nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(&mesh, normalizer);
 
   bc.label(data.label);
   bc.BoundaryConditionPoints::_getPoints(mesh);

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryMesh.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryMesh.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -23,6 +23,7 @@
 #include "data/BoundaryMeshData.hh" // USES BoundaryMeshData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
@@ -76,7 +77,7 @@
   cs.setSpaceDim(mesh.dimension());
   cs.initialize();
   mesh.coordsys(&cs);
-  mesh.nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(&mesh, normalizer);
 
   // Create submesh
   topology::SubMesh submesh(mesh, _data->bcLabel);
@@ -137,7 +138,7 @@
   cs.setSpaceDim(mesh.dimension());
   cs.initialize();
   mesh.coordsys(&cs);
-  mesh.nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(&mesh, normalizer);
 
   // Adjust topology
   faults::FaultCohesiveKin fault;

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -25,6 +25,7 @@
 #include "data/DirichletData.hh" // USES DirichletData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
 #include "pylith/topology/Stratum.hh" // USES Stratum
@@ -525,7 +526,7 @@
   cs.setSpaceDim(mesh->dimension());
   cs.initialize();
   mesh->coordsys(&cs);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   spatialdata::spatialdb::SimpleDB db("TestDirichletBC initial");
   spatialdata::spatialdb::SimpleIOAscii dbIO;

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -25,6 +25,7 @@
 #include "data/DirichletDataMulti.hh" // USES DirichletData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VisitorMesh
@@ -362,7 +363,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Setup boundary condition A
   spatialdata::spatialdb::SimpleDB db("TestDirichletBCMulti initial");

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBoundary.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBoundary.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBoundary.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -25,6 +25,7 @@
 #include "data/DirichletData.hh" // USES DirichletData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
@@ -121,7 +122,7 @@
   cs.setSpaceDim(mesh->dimension());
   cs.initialize();
   mesh->coordsys(&cs);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   spatialdata::spatialdb::SimpleDB db("TestDirichletBoundary initial");
   spatialdata::spatialdb::SimpleIOAscii dbIO;

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -25,6 +25,7 @@
 #include "data/NeumannData.hh" // USES NeumannData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Fields.hh" // USES Fields
@@ -698,7 +699,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Set up quadrature
   _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -25,6 +25,7 @@
 #include "data/PointForceData.hh" // USES PointForceData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
 #include "pylith/topology/SolutionFields.hh" // USES SolutionFields
@@ -259,7 +260,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   spatialdata::spatialdb::SimpleDB dbInitial("TestPointForce initial");
   spatialdata::spatialdb::SimpleIOAscii dbInitialIO;

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependentPoints.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependentPoints.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependentPoints.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -23,6 +23,7 @@
 #include "pylith/bc/PointForce.hh" // USES PointForce
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitormesh
@@ -135,7 +136,7 @@
   normalizer.pressureScale(_TestTimeDependentPoints::pressureScale);
   normalizer.lengthScale(_TestTimeDependentPoints::lengthScale);
   normalizer.timeScale(_TestTimeDependentPoints::timeScale);
-  _mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(_mesh, normalizer);
 
   _bc = new PointForce();CPPUNIT_ASSERT(_bc);
   _bc->label("bc");

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -27,6 +27,7 @@
 #include "pylith/faults/BruneSlipFn.hh" // USES BruneSlipFn
 #include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Stratum.hh" // USES Stratum
@@ -207,7 +208,7 @@
   normalizer.pressureScale(_TestEqKinSrc::pressureScale);
   normalizer.densityScale(_TestEqKinSrc::densityScale);
   normalizer.timeScale(_TestEqKinSrc::timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Create fault mesh
   TestFaultMesh::createFaultMesh(faultMesh, mesh, faultLabel, faultId);

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -26,6 +26,7 @@
 #include "pylith/faults/FaultCohesiveTract.hh" // USES FaultsCohesiveTract
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
 
@@ -691,7 +692,7 @@
   mesh.coordsys(&cs);
 
   spatialdata::units::Nondimensional normalizer;
-  mesh.nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(&mesh, normalizer);
 
   PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
   PetscInt firstFaultVertex = 0;

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -26,6 +26,7 @@
 #include "data/CohesiveDynData.hh" // USES CohesiveDynData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/topology/SolutionFields.hh" // USES SolutionFields
@@ -635,7 +636,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
   
   _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
 			  _data->basisDeriv,

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -25,6 +25,7 @@
 #include "data/CohesiveImpulsesData.hh" // USES CohesiveImpulsesData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
@@ -329,7 +330,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
   
   _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
 			  _data->basisDeriv,

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -27,6 +27,7 @@
 #include "pylith/faults/EqKinSrc.hh" // USES EqKinSrc
 #include "pylith/faults/BruneSlipFn.hh" // USES BruneSlipFn
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
@@ -772,13 +773,11 @@
   // Set scales
   // Most test data is insensitive to the scales because we set the fields directly.
   spatialdata::units::Nondimensional normalizer;
-#if 1 // :TODO: USE SCALES
   normalizer.lengthScale(_data->lengthScale);
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-#endif
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
   
   //mesh->debug(true); // DEBUGGING
   

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -26,6 +26,7 @@
 
 #include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
@@ -348,7 +349,7 @@
   normalizer.pressureScale(_TestTractPerturbation::pressureScale);
   normalizer.densityScale(_TestTractPerturbation::densityScale);
   normalizer.timeScale(_TestTractPerturbation::timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Create fault mesh
   TestFaultMesh::createFaultMesh(faultMesh, mesh, faultLabel, faultId);

Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -26,6 +26,7 @@
 #include "pylith/materials/ElasticIsotropic3D.hh" // USES ElasticIsotropic3D
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
@@ -371,7 +372,7 @@
   const PylithScalar dt = _data->dt;
 
   // Setup mesh
-  mesh->createDMMesh(_data->cellDim);
+  topology::MeshOps::createDMMesh(mesh, _data->cellDim);
   PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
 
   // Cells and vertices
@@ -411,7 +412,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Setup material
   spatialdata::spatialdb::SimpleIOAscii iohandler;

Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitLgDeform.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitLgDeform.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -26,6 +26,7 @@
 #include "pylith/materials/ElasticIsotropic3D.hh" // USES ElasticIsotropic3D
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
@@ -264,7 +265,7 @@
   const PylithScalar dt = _data->dt;
 
   // Setup mesh
-  mesh->createDMMesh(_data->cellDim);
+  topology::MeshOps::createDMMesh(mesh, _data->cellDim);
   PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
 
   // Cells and vertices
@@ -304,7 +305,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Setup material
   spatialdata::spatialdb::SimpleIOAscii iohandler;

Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -27,6 +27,7 @@
 #include "pylith/materials/ElasticIsotropic3D.hh" // USES ElasticIsotropic3D
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
@@ -359,7 +360,7 @@
   const PylithScalar dt = _data->dt;
 
   // Setup mesh
-  mesh->createDMMesh(_data->cellDim);
+  topology::MeshOps::createDMMesh(mesh, _data->cellDim);
   PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
 
   // Cells and vertices
@@ -399,7 +400,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Setup material
   spatialdata::spatialdb::SimpleIOAscii iohandler;

Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -27,6 +27,7 @@
 #include "pylith/materials/ElasticPlaneStrain.hh" // USES ElasticPlaneStrain
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
@@ -358,7 +359,7 @@
   const PylithScalar dt = _data->dt;
 
   // Setup mesh
-  mesh->createDMMesh(_data->cellDim);
+  topology::MeshOps::createDMMesh(mesh, _data->cellDim);
   PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
 
   // Cells and vertices
@@ -398,7 +399,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Setup material
   spatialdata::spatialdb::SimpleIOAscii iohandler;

Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -27,6 +27,7 @@
 #include "pylith/materials/ElasticIsotropic3D.hh" // USES ElasticIsotropic3D
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
@@ -342,7 +343,7 @@
   const int spaceDim = _data->spaceDim;
 
   // Setup mesh
-  mesh->createDMMesh(_data->cellDim);
+  topology::MeshOps::createDMMesh(mesh, _data->cellDim);
   PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
 
   // Cells and vertices
@@ -382,7 +383,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Setup material
   spatialdata::spatialdb::SimpleIOAscii iohandler;

Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicitLgDeform.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicitLgDeform.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -26,6 +26,7 @@
 #include "pylith/materials/ElasticIsotropic3D.hh" // USES ElasticIsotropic3D
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
@@ -264,7 +265,7 @@
   const int spaceDim = _data->spaceDim;
 
   // Setup mesh
-  mesh->createDMMesh(_data->cellDim);
+  topology::MeshOps::createDMMesh(mesh, _data->cellDim);
   PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
 
   // Cells and vertices
@@ -303,7 +304,7 @@
   normalizer.pressureScale(_data->pressureScale);
   normalizer.densityScale(_data->densityScale);
   normalizer.timeScale(_data->timeScale);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Setup material
   spatialdata::spatialdb::SimpleIOAscii iohandler;

Modified: short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -27,6 +27,7 @@
 #include "pylith/friction/SlipWeakening.hh" // USES SlipWeakening
 
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
 #include "pylith/topology/Stratum.hh" // USES Stratum
@@ -818,17 +819,19 @@
   iohandler.filename("data/tri3.mesh");
   iohandler.read(mesh);
 
-  // Set up coordinates
+  // Setup coordinates
   spatialdata::geocoords::CSCart cs;
+  cs.setSpaceDim(mesh->dimension());
+  cs.initialize();
+  mesh->coordsys(&cs);
+  
+  // Setup scales
   spatialdata::units::Nondimensional normalizer;
   normalizer.lengthScale(data->lengthScale);
   normalizer.pressureScale(data->pressureScale);
   normalizer.timeScale(data->timeScale);
   normalizer.densityScale(data->densityScale);
-  cs.setSpaceDim(mesh->dimension());
-  cs.initialize();
-  mesh->coordsys(&cs);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Setup quadrature
   feassemble::Quadrature<topology::SubMesh> quadrature;

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -24,6 +24,7 @@
 #include "data/ElasticStrain1DData.hh" // USES ElasticStrain1DData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Fields.hh" // USES Fields
 #include "pylith/topology/Stratum.hh" // USES StratumIS
@@ -828,17 +829,19 @@
   iohandler.filename("data/line3.mesh");
   iohandler.read(mesh);
 
-  // Set up coordinates
+  // Setup coordinates.
   spatialdata::geocoords::CSCart cs;
+  cs.setSpaceDim(mesh->dimension());
+  cs.initialize();
+  mesh->coordsys(&cs);
+
+  // Setup scales.
   spatialdata::units::Nondimensional normalizer;
   normalizer.lengthScale(data->lengthScale);
   normalizer.pressureScale(data->pressureScale);
   normalizer.timeScale(data->timeScale);
   normalizer.densityScale(data->densityScale);
-  cs.setSpaceDim(mesh->dimension());
-  cs.initialize();
-  mesh->coordsys(&cs);
-  mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(mesh, normalizer);
 
   // Setup quadrature
   feassemble::Quadrature<topology::Mesh> quadrature;

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -25,6 +25,7 @@
 #include "pylith/utils/array.hh" // USES scalar_array
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Stratum.hh" // USES StratumIS
 #include "pylith/topology/VisitorMesh.hh" // USES VisitorMesh
@@ -221,7 +222,7 @@
   normalizer.pressureScale(pressureScale);
   normalizer.timeScale(timeScale);
   normalizer.densityScale(densityScale);
-  mesh.nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(&mesh, normalizer);
 
   // Setup quadrature
   feassemble::Quadrature<topology::Mesh> quadrature;

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterFaultMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterFaultMesh.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterFaultMesh.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -101,7 +101,7 @@
   fault.label(_data->faultLabel);
   fault.id(_data->faultId);
   fault.adjustTopology(_mesh, &firstFaultVertex, &firstLagrangeVertex, &firstFaultCell, _flipFault);
-  faults::CohesiveTopology::createFaultParallel(_faultMesh, *_mesh, _data->faultId, useLagrangeConstraints);
+  faults::CohesiveTopology::createFaultParallel(_faultMesh, *_mesh, _data->faultId, _data->faultLabel, useLagrangeConstraints);
 
   PYLITH_METHOD_END;
 } // _initialize

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -21,6 +21,7 @@
 #include "TestMeshIO.hh" // Implementation of class methods
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::nondimensionalize()
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
 #include "pylith/meshio/MeshIO.hh" // USES MeshIO
@@ -88,7 +89,7 @@
   PetscErrorCode err;
 
   err = DMPlexCreateFromCellList(_mesh->comm(), data.cellDim, data.numCells, data.numVertices, data.numCorners, interpolateMesh, data.cells, data.spaceDim, data.vertices, &dmMesh);PYLITH_CHECK_ERROR(err);
-  _mesh->setDMMesh(dmMesh);
+  _mesh->dmMesh(dmMesh);
 
   // Material ids
   PetscInt cStart, cEnd;
@@ -125,7 +126,7 @@
   cs.setSpaceDim(data.spaceDim);
   cs.initialize();
   _mesh->coordsys(&cs);
-  _mesh->nondimensionalize(normalizer);
+  topology::MeshOps::nondimensionalize(_mesh, normalizer);
 
   PYLITH_METHOD_END;
 } // _createMesh

Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -22,6 +22,7 @@
 
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::createDMMesh()
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
 
@@ -1434,7 +1435,7 @@
 
   PetscErrorCode err = 0;
 
-  mesh->createDMMesh(_TestFieldMesh::cellDim);
+  MeshOps::createDMMesh(mesh, _TestFieldMesh::cellDim);
   PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
   
   err = DMPlexSetChart(dmMesh, 0, ncells+nvertices);PYLITH_CHECK_ERROR(err);

Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldSubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldSubMesh.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldSubMesh.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -22,6 +22,7 @@
 
 #include "pylith/topology/Field.hh" // USES Field
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::createDMMesh()
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
 
@@ -1007,7 +1008,7 @@
 
   PetscErrorCode err = 0;
 
-  mesh->createDMMesh(_TestFieldSubMesh::cellDim);
+  MeshOps::createDMMesh(mesh, _TestFieldSubMesh::cellDim);
   PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);  
   err = DMPlexSetChart(dmMesh, 0, ncells+nvertices);PYLITH_CHECK_ERROR(err);
   for(PetscInt c = 0; c < ncells; ++c) {

Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -19,6 +19,7 @@
 #include <portinfo>
 
 #include "TestMesh.hh" // Implementation of class methods
+
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
 #include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
@@ -26,7 +27,6 @@
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 
 #include "spatialdata/geocoords/CSCart.hh" // USES CSCart
-#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
 // ----------------------------------------------------------------------
 CPPUNIT_TEST_SUITE_REGISTRATION( pylith::topology::TestMesh );
@@ -41,7 +41,7 @@
   int result = 0;
 
   Mesh mesh;
-  CPPUNIT_ASSERT(!mesh._newMesh);
+  CPPUNIT_ASSERT(!mesh._dmMesh);
   CPPUNIT_ASSERT_EQUAL(0, mesh.dimension());
   CPPUNIT_ASSERT_EQUAL(false, mesh.debug());
   MPI_Comm_compare(PETSC_COMM_WORLD, mesh.comm(), &result);
@@ -49,14 +49,14 @@
 
   int dim = 2;
   Mesh mesh2(dim);
-  CPPUNIT_ASSERT(mesh2._newMesh);
+  CPPUNIT_ASSERT(mesh2._dmMesh);
   CPPUNIT_ASSERT_EQUAL(dim, mesh2.dimension());
   MPI_Comm_compare(PETSC_COMM_WORLD, mesh2.comm(), &result);
   CPPUNIT_ASSERT_EQUAL(int(MPI_IDENT), result);
 
   dim = 1;
   Mesh mesh3(dim, PETSC_COMM_SELF);
-  CPPUNIT_ASSERT(mesh3._newMesh);
+  CPPUNIT_ASSERT(mesh3._dmMesh);
   CPPUNIT_ASSERT_EQUAL(dim, mesh3.dimension());
   MPI_Comm_compare(PETSC_COMM_WORLD, mesh3.comm(), &result);
   CPPUNIT_ASSERT_EQUAL(int(MPI_CONGRUENT), result);
@@ -65,29 +65,6 @@
 } // testConstructor
 
 // ----------------------------------------------------------------------
-// Test createDMMesh().
-void
-pylith::topology::TestMesh::testCreateDMMesh(void)
-{ // testCreateDMMesh
-  PYLITH_METHOD_BEGIN;
-
-  Mesh mesh;
-  CPPUNIT_ASSERT(!mesh._newMesh);
-
-  int dim = 2;
-  mesh.createDMMesh(dim);
-  CPPUNIT_ASSERT(mesh._newMesh);
-  CPPUNIT_ASSERT_EQUAL(dim, mesh.dimension());
-
-  dim = 1;
-  mesh.createDMMesh(dim);
-  CPPUNIT_ASSERT(mesh._newMesh);
-  CPPUNIT_ASSERT_EQUAL(dim, mesh.dimension());
-
-  PYLITH_METHOD_END;
-} // testCreateDMMesh
-
-// ----------------------------------------------------------------------
 // Test dmMesh().
 void
 pylith::topology::TestMesh::testDMMesh(void)
@@ -98,8 +75,7 @@
   PetscInt dmDim;
   Mesh mesh(dim);
   
-  DM dmMesh = mesh.dmMesh();
-  CPPUNIT_ASSERT(dmMesh);
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
   PetscErrorCode err = DMPlexGetDimension(dmMesh, &dmDim);PYLITH_CHECK_ERROR(err);
   CPPUNIT_ASSERT_EQUAL(dim, dmDim);
 
@@ -174,61 +150,4 @@
   PYLITH_METHOD_END;
 } // testComm
 
-// ----------------------------------------------------------------------
-// Test nondimensionalize().
-void
-pylith::topology::TestMesh::testNondimensionalize(void)
-{ // testNondimensionalizer
-  PYLITH_METHOD_BEGIN;
-
-  const PylithScalar lengthScale = 2.0;
-  const int spaceDim = 2;
-  const int numVertices = 4;
-  const int coordinates[] = { 
-    -1.0, 0.0,
-    0.0, -1.0,
-    0.0, 1.0,
-    1.0, 0.0,
-  };
-
-  Mesh mesh;
-  meshio::MeshIOAscii iohandler;
-  iohandler.filename("data/tri3.mesh");
-  iohandler.read(&mesh);
-
-  spatialdata::geocoords::CSCart cs;
-  cs.setSpaceDim(2);
-  mesh.coordsys(&cs);
-  spatialdata::units::Nondimensional normalizer;
-  normalizer.lengthScale(lengthScale);
-  mesh.nondimensionalize(normalizer);
-
-  // Get vertices
-  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
-  Stratum depthStratum(dmMesh, Stratum::DEPTH, 0);
-  const PetscInt vStart = depthStratum.begin();
-  const PetscInt vEnd = depthStratum.end();
-
-  // Check nondimensional coordinates
-  CoordsVisitor coordsVisitor(dmMesh);
-  const PetscScalar* coordsArray = coordsVisitor.localArray();
-
-  const PylithScalar tolerance = 1.0e-06;
-  for(PetscInt v = vStart, i = 0; v < vEnd; ++v) {
-    CPPUNIT_ASSERT_EQUAL(spaceDim, coordsVisitor.sectionDof(v));
-    const PetscInt off = coordsVisitor.sectionOffset(v);
-    for (int iDim=0; iDim < spaceDim; ++iDim, ++i) {
-      const PylithScalar coordE = coordinates[i] / lengthScale;
-      if (fabs(coordE) < 1.0) {
-        CPPUNIT_ASSERT_DOUBLES_EQUAL(coordE, coordsArray[off+iDim], tolerance);
-      } else {
-        CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, coordsArray[off+iDim]/coordE, tolerance);
-      } // if/else
-    } // for
-  } // for
-  
-  PYLITH_METHOD_END;
-} // testNondimensionalize
-
-
 // End of file 

Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.hh	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.hh	2013-05-13 18:43:34 UTC (rev 22057)
@@ -47,13 +47,11 @@
   CPPUNIT_TEST_SUITE( TestMesh );
 
   CPPUNIT_TEST( testConstructor );
-  CPPUNIT_TEST( testCreateDMMesh );
   CPPUNIT_TEST( testDMMesh );
   CPPUNIT_TEST( testCoordsys );
   CPPUNIT_TEST( testDebug );
   CPPUNIT_TEST( testDimension );
   CPPUNIT_TEST( testComm );
-  CPPUNIT_TEST( testNondimensionalize );
 
   CPPUNIT_TEST_SUITE_END();
 
@@ -63,9 +61,6 @@
   /// Test constructor.
   void testConstructor(void);
 
-  /// Test createDMMesh().
-  void testCreateDMMesh(void);
-
   /// Test dmMesh().
   void testDMMesh(void);
 
@@ -81,9 +76,6 @@
   /// Test comm().
   void testComm(void);
 
-  /// Test nondimensionalize().
-  void testNondimensionalize(void);
-
 }; // class TestMesh
 
 #endif // pylith_topology_testmesh_hh

Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -23,13 +23,100 @@
 #include "pylith/topology/MeshOps.hh" // USES MeshOps
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
 #include <stdexcept> // USES std::runtime_error
 
 // ----------------------------------------------------------------------
 CPPUNIT_TEST_SUITE_REGISTRATION( pylith::topology::TestMeshOps );
 
 // ----------------------------------------------------------------------
+// Test createDMMesh().
+void
+pylith::topology::TestMeshOps::testCreateDMMesh(void)
+{ // testCreateDMMesh
+  PYLITH_METHOD_BEGIN;
+
+  Mesh mesh;
+  CPPUNIT_ASSERT(!mesh.dmMesh());
+
+  int dim = 2;
+  MeshOps::createDMMesh(&mesh, dim, PETSC_COMM_SELF, "domain");
+  CPPUNIT_ASSERT(mesh.dmMesh());
+  CPPUNIT_ASSERT_EQUAL(dim, mesh.dimension());
+
+  dim = 1;
+  MeshOps::createDMMesh(&mesh, dim);
+  CPPUNIT_ASSERT(mesh.dmMesh());
+  CPPUNIT_ASSERT_EQUAL(dim, mesh.dimension());
+
+  PYLITH_METHOD_END;
+} // testCreateDMMesh
+
+
+// ----------------------------------------------------------------------
+// Test nondimensionalize().
+void
+pylith::topology::TestMeshOps::testNondimensionalize(void)
+{ // testNondimensionalize
+  PYLITH_METHOD_BEGIN;
+
+  const PylithScalar lengthScale = 2.0;
+  const int spaceDim = 2;
+  const int numVertices = 4;
+  const int coordinates[numVertices*spaceDim] = { 
+    -1.0, 0.0,
+    0.0, -1.0,
+    0.0, 1.0,
+    1.0, 0.0,
+  };
+
+  Mesh mesh;
+  meshio::MeshIOAscii iohandler;
+  iohandler.filename("data/tri3.mesh");
+  iohandler.read(&mesh);
+
+  spatialdata::geocoords::CSCart cs;
+  cs.setSpaceDim(2);
+  mesh.coordsys(&cs);
+  spatialdata::units::Nondimensional normalizer;
+  normalizer.lengthScale(lengthScale);
+  MeshOps::nondimensionalize(&mesh, normalizer);
+
+  // Get vertices
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  Stratum depthStratum(dmMesh, Stratum::DEPTH, 0);
+  const PetscInt vStart = depthStratum.begin();
+  const PetscInt vEnd = depthStratum.end();
+
+  // Check nondimensional coordinates
+  CoordsVisitor coordsVisitor(dmMesh);
+  const PetscScalar* coordsArray = coordsVisitor.localArray();CPPUNIT_ASSERT(coordsArray);
+
+  const PylithScalar tolerance = 1.0e-06;
+  for(PetscInt v = vStart, i = 0; v < vEnd; ++v) {
+    CPPUNIT_ASSERT_EQUAL(spaceDim, coordsVisitor.sectionDof(v));
+    const PetscInt off = coordsVisitor.sectionOffset(v);
+    for (int iDim=0; iDim < spaceDim; ++iDim, ++i) {
+      const PylithScalar coordE = coordinates[i] / lengthScale;
+      if (fabs(coordE) < 1.0) {
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(coordE, coordsArray[off+iDim], tolerance);
+      } else {
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, coordsArray[off+iDim]/coordE, tolerance);
+      } // if/else
+    } // for
+  } // for
+  
+  PYLITH_METHOD_END;
+} // testNondimensionalize
+
+
+// ----------------------------------------------------------------------
 // Test checkMaterialIds().
 void
 pylith::topology::TestMeshOps::testCheckMaterialIds(void)

Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.hh	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.hh	2013-05-13 18:43:34 UTC (rev 22057)
@@ -43,6 +43,8 @@
   // CPPUNIT TEST SUITE /////////////////////////////////////////////////
   CPPUNIT_TEST_SUITE( TestMeshOps );
 
+  CPPUNIT_TEST( testCreateDMMesh );
+  CPPUNIT_TEST( testNondimensionalize );
   CPPUNIT_TEST( testCheckMaterialIds );
 
   CPPUNIT_TEST_SUITE_END();
@@ -50,6 +52,12 @@
   // PUBLIC METHODS /////////////////////////////////////////////////////
 public :
 
+  /// Test createDMMesh().
+  void testCreateDMMesh(void);
+
+  /// Test nondimensionalize().
+  void testNondimensionalize(void);
+
   /// Test checkMaterialIds().
   void testCheckMaterialIds(void);
 

Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.cc	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.cc	2013-05-13 18:43:34 UTC (rev 22057)
@@ -23,6 +23,7 @@
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/MeshOps.hh" // USES MeshOps::createDMMesh()
 
 #include "spatialdata/geocoords/CSCart.hh" // USES CSCart
 
@@ -284,23 +285,7 @@
 } // testComm
 
 // ----------------------------------------------------------------------
-// Test initialize().
 void
-pylith::topology::TestSubMesh::testInitialize(void)
-{ // testInitialize
-  PYLITH_METHOD_BEGIN;
-
-  Mesh mesh2D;
-  _buildMesh(&mesh2D);
-  SubMesh mesh(mesh2D, _TestSubMesh::label);
-
-  mesh.initialize();
-
-  PYLITH_METHOD_END;
-} // testInitialize
-
-// ----------------------------------------------------------------------
-void
 pylith::topology::TestSubMesh::_buildMesh(Mesh* mesh)
 { // _buildMesh
   PYLITH_METHOD_BEGIN;
@@ -316,7 +301,7 @@
   const PylithScalar* coordinates = _TestSubMesh::coordinates;
   const bool interpolate = false;
 
-  mesh->createDMMesh(_TestSubMesh::cellDim);
+  MeshOps::createDMMesh(mesh, _TestSubMesh::cellDim);
   PetscDM dmMesh = mesh->dmMesh();CPPUNIT_ASSERT(dmMesh);
   PetscErrorCode err;
   

Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.hh	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.hh	2013-05-13 18:43:34 UTC (rev 22057)
@@ -56,7 +56,6 @@
   CPPUNIT_TEST( testNumVertices );
   CPPUNIT_TEST( testNumCells );
   CPPUNIT_TEST( testComm );
-  CPPUNIT_TEST( testInitialize );
 
   CPPUNIT_TEST_SUITE_END();
 
@@ -93,9 +92,6 @@
   /// Test comm().
   void testComm(void);
 
-  /// Test initialize().
-  void testInitialize(void);
-
 // PRIVATE METHODS /////////////////////////////////////////////////////
 private :
 

Modified: short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveDyn.py	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveDyn.py	2013-05-13 18:43:34 UTC (rev 22057)
@@ -315,9 +315,9 @@
     firstFaultCell      = 2*nvertices
     fault.adjustTopology(mesh, firstFaultVertex, firstLagrangeVertex,
                          firstFaultCell)
+    from pylith.topology.topology import MeshOps_nondimensionalize
+    MeshOps_nondimensionalize(mesh, normalizer)
 
-    mesh.nondimensionalize(normalizer)
-
     fault.preinitialize(mesh)
     fault.timeStep(dt)
     fault.verifyConfiguration()

Modified: short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py	2013-05-13 18:43:34 UTC (rev 22057)
@@ -63,21 +63,6 @@
     return
 
 
-  def test_createDMMesh(self):
-    """
-    Test createDMMesh().
-    """
-    mesh = Mesh()
-
-    mesh.createDMMesh()
-    self.assertEqual(3, mesh.dimension())
-
-    dim = 2
-    mesh.createDMMesh(dim)
-    self.assertEqual(dim, mesh.dimension())
-    return
-
-
   def test_coordsys(self):
     """
     Test coordsys().

Modified: short/3D/PyLith/trunk/unittests/pytests/topology/TestSubMesh.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/topology/TestSubMesh.py	2013-05-13 16:03:35 UTC (rev 22056)
+++ short/3D/PyLith/trunk/unittests/pytests/topology/TestSubMesh.py	2013-05-13 18:43:34 UTC (rev 22057)
@@ -110,17 +110,6 @@
     return
 
 
-  def test_initialize(self):
-    """
-    Test initialize().
-    """
-    mesh = self._getMesh()
-    submesh = SubMesh(mesh, "bc")
-
-    submesh.initialize()
-    return
-
-
   # PRIVATE METHODS ////////////////////////////////////////////////////
 
   def _getMesh(self):



More information about the CIG-COMMITS mailing list