[cig-commits] r20803 - in short/3D/PyLith/trunk: libsrc/pylith/faults libsrc/pylith/feassemble libsrc/pylith/materials libsrc/pylith/meshio libsrc/pylith/problems libsrc/pylith/topology unittests/libtests/friction unittests/libtests/meshio

knepley at geodynamics.org knepley at geodynamics.org
Thu Oct 4 05:17:01 PDT 2012


Author: knepley
Date: 2012-10-04 05:17:00 -0700 (Thu, 04 Oct 2012)
New Revision: 20803

Modified:
   short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
   short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIO.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/VertexFilterVecNorm.cc
   short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/FieldBase.hh
   short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc
   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/unittests/libtests/friction/TestFrictionModel.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestCellFilterAvg.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterBCMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterFaultMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterPoints.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterSubMesh.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestOutputManager.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestVertexFilterVecNorm.cc
Log:
Updating tests in meshio, now only points tests have old sections

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -775,15 +775,14 @@
   SieveSubMesh::label_sequence::const_iterator fVertices2End   = fVertices2->end();
 
   PetscSection coordSection, newCoordSection;
-  Vec          coordinatesDM, newCoordinatesDM;
+  Vec          coordinatesVec, newCoordinatesVec;
   PetscScalar *coords, *newCoords;
   PetscInt     coordSize;
  
 #ifdef USE_DMCOMPLEX_ON
   err = DMComplexGetCoordinateSection(complexMesh, &coordSection);CHECK_PETSC_ERROR(err);
   err = DMComplexGetCoordinateSection(newMesh,     &newCoordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(complexMesh, &coordinatesDM);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(newMesh,     &newCoordinatesDM);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(complexMesh, &coordinatesVec);CHECK_PETSC_ERROR(err);
   err = PetscSectionSetChart(newCoordSection, vStart+extraCells, vEnd+extraCells+extraVertices);CHECK_PETSC_ERROR(err);
   for(PetscInt v = vStart; v < vEnd; ++v) {
     PetscInt dof;
@@ -807,10 +806,11 @@
 #ifdef USE_DMCOMPLEX_ON
   err = PetscSectionSetUp(newCoordSection);CHECK_PETSC_ERROR(err);
   err = PetscSectionGetStorageSize(newCoordSection, &coordSize);CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(newCoordinatesDM, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(newCoordinatesDM);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(coordinatesDM, &coords);CHECK_PETSC_ERROR(err);
-  err = VecGetArray(newCoordinatesDM, &newCoords);CHECK_PETSC_ERROR(err);
+  err = VecCreate(((PetscObject) newMesh)->comm, &newCoordinatesVec);CHECK_PETSC_ERROR(err);
+  err = VecSetSizes(newCoordinatesVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+  err = VecSetFromOptions(newCoordinatesVec);CHECK_PETSC_ERROR(err);
+  err = VecGetArray(coordinatesVec, &coords);CHECK_PETSC_ERROR(err);
+  err = VecGetArray(newCoordinatesVec, &newCoords);CHECK_PETSC_ERROR(err);
 
   for(PetscInt v = vStart; v < vEnd; ++v) {
     PetscInt dof, off, newoff, d;
@@ -847,13 +847,13 @@
     } // if
 #endif
   } // for
-  //err = VecRestoreArray(coordinatesDM, &coords);CHECK_PETSC_ERROR(err);
-  //err = VecRestoreArray(newCoordinatesDM, &newCoords);CHECK_PETSC_ERROR(err);
+  err = VecRestoreArray(coordinatesVec, &coords);CHECK_PETSC_ERROR(err);
+  err = VecRestoreArray(newCoordinatesVec, &newCoords);CHECK_PETSC_ERROR(err);
+  err = DMSetCoordinatesLocal(newMesh, newCoordinatesVec);CHECK_PETSC_ERROR(err);
   if (debug) coordinates->view("Coordinates with shadow vertices");
 
   logger.stagePop();
 
-  err = DMDestroy(&complexMesh);CHECK_PETSC_ERROR(err);
   mesh->setDMMesh(newMesh);
 } // create
 
@@ -876,6 +876,8 @@
 
   const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
   assert(!sieveMesh.isNull());
+  DM dmMesh = mesh.dmMesh();
+  assert(dmMesh);
   ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
 
   const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
@@ -967,6 +969,11 @@
 					     fRenumbering,
 					     fault->getArrowSection("orientation").ptr());
   }
+  SieveMesh::renumbering_type convertRenumbering;
+  DM dmFaultMesh;
+
+  ALE::ISieveConverter::convertMesh(*fault, &dmFaultMesh, convertRenumbering);
+  faultMesh->setDMMesh(dmFaultMesh);
   fault      = NULL;
   faultSieve = NULL;
 
@@ -1016,7 +1023,47 @@
 			      coordinates->restrictPoint(*v_iter));
   }
   //faultSieveMesh->view("Parallel fault mesh");
+  const SieveMesh::renumbering_type::const_iterator convertRenumberingEnd = convertRenumbering.end();
+  PetscSection   coordSection, faultCoordSection;
+  Vec            coordinateVec, faultCoordinateVec;
+  PetscScalar   *a, *fa;
+  PetscInt       pvStart, pvEnd, fvStart, fvEnd, n;
+  PetscErrorCode err;
 
+  err = DMComplexGetDepthStratum(dmMesh, 0, &pvStart, &pvEnd);CHECK_PETSC_ERROR(err);
+  err = DMComplexGetDepthStratum(dmFaultMesh, 0, &fvStart, &fvEnd);CHECK_PETSC_ERROR(err);
+  err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+  err = DMComplexGetCoordinateSection(dmFaultMesh, &faultCoordSection);CHECK_PETSC_ERROR(err);
+  err = PetscSectionSetChart(faultCoordSection, fvStart, fvEnd);CHECK_PETSC_ERROR(err);
+  for(PetscInt v = pvStart; v < pvEnd; ++v) {
+    PetscInt dof;
+
+    if (convertRenumbering.find(v) == convertRenumberingEnd) continue;
+    err = PetscSectionGetDof(coordSection, v, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionSetDof(faultCoordSection, convertRenumbering[v], dof);CHECK_PETSC_ERROR(err);
+  }
+  err = PetscSectionSetUp(faultCoordSection);CHECK_PETSC_ERROR(err);
+  err = PetscSectionGetStorageSize(faultCoordSection, &n);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordinateVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmFaultMesh, &faultCoordinateVec);CHECK_PETSC_ERROR(err);
+  err = VecSetSizes(faultCoordinateVec, n, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+  err = VecSetFromOptions(faultCoordinateVec);CHECK_PETSC_ERROR(err);
+  err = VecGetArray(coordinateVec, &a);CHECK_PETSC_ERROR(err);
+  err = VecGetArray(faultCoordinateVec, &fa);CHECK_PETSC_ERROR(err);
+  for(PetscInt v = pvStart; v < pvEnd; ++v) {
+    PetscInt dof, off, foff;
+
+    if (convertRenumbering.find(v) == convertRenumberingEnd) continue;
+    err = PetscSectionGetDof(coordSection, v, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(coordSection, v, &off);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(faultCoordSection, convertRenumbering[v], &foff);CHECK_PETSC_ERROR(err);
+    for(PetscInt d = 0; d < dof; ++d) {
+      fa[foff+d] = a[off+d];
+    }
+  }
+  err = VecRestoreArray(coordinateVec, &a);CHECK_PETSC_ERROR(err);
+  err = VecRestoreArray(faultCoordinateVec, &fa);CHECK_PETSC_ERROR(err);
+
   // Update dimensioned coordinates if they exist.
   if (sieveMesh->hasRealSection("coordinates_dimensioned")) {
     const ALE::Obj<topology::Mesh::RealSection>& coordinatesDim =

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -1316,16 +1316,17 @@
   topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
   const topology::Field<topology::SubMesh>& dispRel = 
     _fields->get("relative disp");
+  orientation.addField("orientation", cohesiveDim+1);
+  orientation.setupFields();
   orientation.newSection(dispRel, orientationSize);
-  const ALE::Obj<RealSection>& orientationSection = orientation.section();
-  assert(!orientationSection.isNull());
-  // Create subspaces for along-strike, up-dip, and normal directions
-  for (int iDim = 0; iDim <= cohesiveDim; ++iDim)
-    orientationSection->addSpace();
-  for (int iDim = 0; iDim <= cohesiveDim; ++iDim)
-    orientationSection->setFiberDimension(vertices, spaceDim, iDim);
+  // Create components for along-strike, up-dip, and normal directions
+  orientation.updateDof("orientation", pylith::topology::FieldBase::POINTS_FIELD, spaceDim);
   orientation.allocate();
   orientation.zero();
+  PetscSection orientationSection = orientation.petscSection();
+  Vec          orientationVec     = orientation.localVector();
+  PetscScalar *orientationArray;
+  PetscErrorCode err;
 
   logger.stagePop();
 
@@ -1362,6 +1363,7 @@
   ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type>
     ncV(*sieve, closureSize);
 
+  err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
   for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin; c_iter
       != cellsEnd; ++c_iter) {
     // Get orientations at fault cell's vertices.
@@ -1386,9 +1388,15 @@
         up);
 
       // Update orientation
-      orientationSection->updateAddPoint(cone[v], &orientationVertex[0]);
+      PetscInt off;
+
+      err = PetscSectionGetOffset(orientationSection, cone[v], &off);CHECK_PETSC_ERROR(err);
+      for(PetscInt d = 0; d < orientationSize; ++d) {
+        orientationArray[off+d] += orientationVertex[d];
+      }
     } // for
   } // for
+  err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
 
   //orientation.view("ORIENTATION BEFORE COMPLETE");
 
@@ -1398,11 +1406,15 @@
   // Loop over vertices, make orientation information unit magnitude
   scalar_array vertexDir(orientationSize);
   int count = 0;
+  err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
   for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
       != verticesEnd; ++v_iter, ++count) {
-    orientationVertex = 0.0;
-    orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
-      orientationVertex.size());
+    PetscInt off;
+
+    err = PetscSectionGetOffset(orientationSection, *v_iter, &off);CHECK_PETSC_ERROR(err);
+    for(PetscInt d = 0; d < orientationSize; ++d) {
+      orientationVertex[d] = orientationArray[off+d];
+    }
     for (int iDim = 0; iDim < spaceDim; ++iDim) {
       PylithScalar mag = 0;
       for (int jDim = 0, index = iDim * spaceDim; jDim < spaceDim; ++jDim)
@@ -1412,11 +1424,14 @@
       for (int jDim = 0, index = iDim * spaceDim; jDim < spaceDim; ++jDim)
         orientationVertex[index + jDim] /= mag;
     } // for
-
-    orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
+    for(PetscInt d = 0; d < orientationSize; ++d) {
+      orientationArray[off+d] = orientationVertex[d];
+    }
   } // for
+  err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
   PetscLogFlops(count * orientationSize * 4);
 
+  err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
   if (1 == cohesiveDim && vertices->size() > 0) {
     // Default sense of positive slip is left-lateral and
     // fault-opening.
@@ -1438,9 +1453,13 @@
     // normal direction because it is tied to how the cohesive cells
     // are created.
     assert(vertices->size() > 0);
-    orientationSection->restrictPoint(*vertices->begin(),
-      &orientationVertex[0], orientationVertex.size());
+    PetscInt off;
 
+    err = PetscSectionGetOffset(orientationSection, *vertices->begin(), &off);CHECK_PETSC_ERROR(err);
+    for(PetscInt d = 0; d < orientationSize; ++d) {
+      orientationVertex[d] = orientationArray[off+d];
+    }
+
     assert(2 == spaceDim);
     const PylithScalar* shearDirVertex = &orientationVertex[0];
     const PylithScalar* normalDirVertex = &orientationVertex[2];
@@ -1456,14 +1475,21 @@
       for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin; 
 	   v_iter != verticesEnd;
 	   ++v_iter) {
-        orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
-					  orientationVertex.size());
-        assert(4 == orientationSection->getFiberDimension(*v_iter));
+        PetscInt dof;
+
+        err = PetscSectionGetDof(orientationSection, *v_iter, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(orientationSection, *v_iter, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < orientationSize; ++d) {
+          orientationVertex[d] = orientationArray[off+d];
+        }
+        assert(4 == dof);
         for (int iDim = 0; iDim < 2; ++iDim) // flip shear
           orientationVertex[ishear + iDim] *= -1.0;
 	
         // Update orientation
-        orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
+        for(PetscInt d = 0; d < orientationSize; ++d) {
+          orientationArray[off+d] = orientationVertex[d];
+        }
       } // for
       PetscLogFlops(3 + count * 2);
     } // if
@@ -1486,11 +1512,14 @@
     // are used are the opposite of what we would want, but we cannot
     // flip the fault normal direction because it is tied to how the
     // cohesive cells are created.
-
     assert(vertices->size() > 0);
-    orientationSection->restrictPoint(*vertices->begin(),
-      &orientationVertex[0], orientationVertex.size());
+    PetscInt off;
 
+    err = PetscSectionGetOffset(orientationSection, *vertices->begin(), &off);CHECK_PETSC_ERROR(err);
+    for(PetscInt d = 0; d < orientationSize; ++d) {
+      orientationVertex[d] = orientationArray[off+d];
+    }
+
     assert(3 == spaceDim);
     const PylithScalar* dipDirVertex = &orientationVertex[3];
     const PylithScalar* normalDirVertex = &orientationVertex[6];
@@ -1515,18 +1544,26 @@
       // Flip dip direction
       for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
 	     != verticesEnd; ++v_iter) {
-        orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
-					  orientationVertex.size());
-        assert(9 == orientationSection->getFiberDimension(*v_iter));
+        PetscInt dof;
+
+        err = PetscSectionGetDof(orientationSection, *v_iter, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(orientationSection, *v_iter, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < orientationSize; ++d) {
+          orientationVertex[d] = orientationArray[off+d];
+        }
+        assert(9 == dof);
         for (int iDim = 0; iDim < 3; ++iDim) // flip dip
           orientationVertex[idip + iDim] *= -1.0;
 	
         // Update direction
-        orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
+        for(PetscInt d = 0; d < orientationSize; ++d) {
+          orientationArray[off+d] = orientationVertex[d];
+        }
       } // for
       PetscLogFlops(5 + count * 3);
     } // if
   } // if
+  err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
 
   //orientation.view("ORIENTATION");
 } // _calcOrientation

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -226,7 +226,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   assert(0 != _normalizer);
@@ -499,7 +499,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   assert(0 != _normalizer);
@@ -724,7 +724,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   // Get sparse matrix
@@ -884,7 +884,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   _logger->eventEnd(setupEvent);

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -221,7 +221,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   assert(0 != _normalizer);
@@ -469,7 +469,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   assert(0 != _normalizer);
@@ -658,7 +658,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   // Get sparse matrix
@@ -792,7 +792,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   _logger->eventEnd(setupEvent);

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -204,7 +204,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   assert(0 != _normalizer);
@@ -515,7 +515,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   assert(0 != _normalizer);
@@ -814,7 +814,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   _logger->eventEnd(setupEvent);
@@ -954,7 +954,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   _logger->eventEnd(setupEvent);

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -204,7 +204,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   assert(0 != _normalizer);
@@ -480,7 +480,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   assert(0 != _normalizer);
@@ -729,7 +729,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   _logger->eventEnd(setupEvent);
@@ -869,7 +869,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   _logger->eventEnd(setupEvent);

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -204,7 +204,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 #endif
 
@@ -406,7 +406,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   _logger->eventEnd(setupEvent);

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -200,7 +200,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   assert(0 != _normalizer);
@@ -414,7 +414,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 
   _logger->eventEnd(setupEvent);

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -215,7 +215,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 #endif
 
@@ -563,7 +563,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 #endif
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -375,7 +375,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 #endif
 
@@ -528,7 +528,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 #endif
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -156,7 +156,7 @@
   PetscSection coordSection;
   Vec          coordVec;
   err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
   assert(coordSection);assert(coordVec);
 #endif
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -123,7 +123,7 @@
       const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
       assert(cs);
       err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-      err = DMComplexGetCoordinateVec(dmMesh, &coordinates);CHECK_PETSC_ERROR(err);
+      err = DMGetCoordinatesLocal(dmMesh, &coordinates);CHECK_PETSC_ERROR(err);
       err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
       err = DMComplexGetVTKBounds(dmMesh, PETSC_NULL, &vMax);CHECK_PETSC_ERROR(err);
       if (vMax >= 0) {vEnd = PetscMin(vEnd, vMax);}

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -92,7 +92,9 @@
     DataWriter<mesh_type, field_type>::open(mesh, numTimeSteps, label, labelId);
     const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
 
+#if 0
     const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
+#endif
     DM dmMesh = mesh.dmMesh();
     assert(dmMesh);
     PetscMPIInt    commRank;
@@ -113,13 +115,14 @@
       H5T_IEEE_F64BE : H5T_IEEE_F32BE;
 
     // Write vertex coordinates
+    const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+    assert(cs);
+#if 0
     const ALE::Obj<typename mesh_type::RealSection>& coordinatesSection = 
       sieveMesh->hasRealSection("coordinates_dimensioned") ?
       sieveMesh->getRealSection("coordinates_dimensioned") :
       sieveMesh->getRealSection("coordinates");
     assert(!coordinatesSection.isNull());
-    const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
-    assert(cs);
     topology::FieldBase::Metadata metadata;
     // :KLUDGE: We would like to use field_type for the coordinates
     // field. However, the mesh coordinates are Field<mesh_type> and
@@ -136,7 +139,20 @@
     coordinates.scatterSectionToVector(context);
     PetscVec coordinatesVector = coordinates.vector(context);
     assert(coordinatesVector);
+#else
+    DM       coordDM;
+    Vec      coordVec, coordinatesVector;
+    PetscInt globalSize;
 
+    /* Should use the coordinate DM here (once I put it in) */
+    err = DMGetCoordinateDM(dmMesh, &coordDM);CHECK_PETSC_ERROR(err);
+    err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+    err = DMGetGlobalVector(coordDM, &coordinatesVector);CHECK_PETSC_ERROR(err);
+    err = VecGetSize(coordinatesVector, &globalSize);CHECK_PETSC_ERROR(err);
+    err = DMLocalToGlobalBegin(coordDM, coordVec, INSERT_VALUES, coordinatesVector);CHECK_PETSC_ERROR(err);
+    err = DMLocalToGlobalEnd(coordDM, coordVec, INSERT_VALUES, coordinatesVector);CHECK_PETSC_ERROR(err);
+#endif
+
     const std::string& filenameVertices = _datasetFilename("vertices");
     err = PetscViewerBinaryOpen(((PetscObject) dmMesh)->comm, filenameVertices.c_str(),
 				FILE_MODE_WRITE,
@@ -146,12 +162,14 @@
     CHECK_PETSC_ERROR(err);
     err = VecView(coordinatesVector, binaryViewer); CHECK_PETSC_ERROR(err);
     err = PetscViewerDestroy(&binaryViewer); CHECK_PETSC_ERROR(err);
+
+    err = DMRestoreGlobalVector(coordDM, &coordinatesVector);CHECK_PETSC_ERROR(err);
     
     // Create external dataset for coordinates    
     if (!commRank) {
       const hsize_t ndims = 2;
       hsize_t dims[ndims];
-      dims[0] = vNumbering->getGlobalSize();
+      dims[0] = globalSize;
       dims[1] = cs->spaceDim();
       _h5->createDatasetRawExternal("/geometry", "vertices", 
 				    filenameVertices.c_str(),
@@ -161,6 +179,7 @@
     // Write cells
 
     // Account for censored cells
+#if 0
     int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
     int cellDepth = 0;
     err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX, 
@@ -199,56 +218,103 @@
     int k = 0;
     const typename label_sequence::const_iterator cellsEnd = cells->end();
     for (typename label_sequence::iterator c_iter=cells->begin();
-	 c_iter != cellsEnd;
-	 ++c_iter)
+         c_iter != cellsEnd;
+         ++c_iter) {
       if (cNumbering->isLocal(*c_iter)) {
-	ncV.clear();
-	ALE::ISieveTraversal<sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
-	const typename ALE::ISieveVisitor::NConeRetriever<sieve_type>::oriented_point_type* cone =
-	  ncV.getOrientedPoints();
-	const int coneSize = ncV.getOrientedSize();
-	if (coneSize != numCorners) {
-	  std::ostringstream msg;
-	  msg << "Inconsistency in topology found for mesh '"
-	      << sieveMesh->getName() << "' during output.\n"
-	      << "Number of vertices (" << coneSize << ") in cell '"
-	      << *c_iter << "' does not expected number of vertices ("
-	      << numCorners << ").";
-	  throw std::runtime_error(msg.str());
-	} // if
-	for (int c=0; c < coneSize; ++c)
-	  tmpVertices[k++] = vNumbering->getIndex(cone[c].first);
+        ncV.clear();
+        ALE::ISieveTraversal<sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
+        const typename ALE::ISieveVisitor::NConeRetriever<sieve_type>::oriented_point_type* cone =
+          ncV.getOrientedPoints();
+        const int coneSize = ncV.getOrientedSize();
+        if (coneSize != numCorners) {
+          const char *name;
+          std::ostringstream msg;
+          err = PetscObjectGetName((PetscObject) dmMesh, &name);CHECK_PETSC_ERROR(err);
+          msg << "Inconsistency in topology found for mesh '" << name << "' during output.\n"
+              << "Number of vertices (" << coneSize << ") in cell '" << c
+              << "' does not expected number of vertices (" << numCorners << ").";
+          throw std::runtime_error(msg.str());
+        } // if
+        for(PetscInt c=0; c < coneSize; ++c)
+          tmpVertices[k++] = vNumbering->getIndex(cone[c].first);
       } // if
+    }
+#else
+    const PetscInt *cells, *vertices;
+    IS              globalCellNumbers, globalVertexNumbers;
+    PetscInt        numCornersLocal = 0;
+    PetscInt        numCorners      = 0;
+    PetscInt        numCellsLocal   = 0, numCells;
+    PetscInt        cStart, cEnd, cMax, vStart, numIndices;
 
+    err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+    err = DMComplexGetHeightStratum(dmMesh, 0, &vStart, PETSC_NULL);CHECK_PETSC_ERROR(err);
+    err = DMComplexGetVTKBounds(dmMesh, &cMax, PETSC_NULL);CHECK_PETSC_ERROR(err);
+    if (cMax >= 0) {cEnd = PetscMin(cEnd, cMax);}
+
+    err = DMComplexGetCellNumbering(dmMesh, &globalCellNumbers);CHECK_PETSC_ERROR(err);
+    err = DMComplexGetVertexNumbering(dmMesh, &globalVertexNumbers);CHECK_PETSC_ERROR(err);
+    if (cEnd-cStart) {err = DMComplexGetConeSize(dmMesh, cStart, &numCornersLocal);CHECK_PETSC_ERROR(err);}
+    err = MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
+    err = MPI_Allreduce(&numCellsLocal,   &numCells,   1, MPIU_INT, MPI_SUM, ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
+    err = ISGetSize(globalCellNumbers, &numIndices);CHECK_PETSC_ERROR(err);
+    err = ISGetIndices(globalCellNumbers, &cells);CHECK_PETSC_ERROR(err);
+    err = ISGetIndices(globalVertexNumbers, &vertices);CHECK_PETSC_ERROR(err);
+    for(PetscInt i = 0; i < numIndices; ++i) {
+      if (cells[i] >= 0) ++numCellsLocal;
+    }
+
+    PetscScalar *tmpVertices = 0;
+    PetscInt     conesSize   = numCellsLocal*numCorners;
+
+    err = PetscMalloc(conesSize * sizeof(PetscScalar), &tmpVertices);CHECK_PETSC_ERROR(err);
+    for(PetscInt c = cStart, k = 0; c < cEnd; ++c) {
+      const PetscInt *cone;
+      PetscInt        coneSize;
+
+      if (cells[c-cStart] < 0) continue;
+      err = DMComplexGetCone(dmMesh, c, &cone);CHECK_PETSC_ERROR(err);
+      err = DMComplexGetConeSize(dmMesh, c, &coneSize);CHECK_PETSC_ERROR(err);
+      if (coneSize != numCorners) {
+        const char *name;
+        std::ostringstream msg;
+        err = PetscObjectGetName((PetscObject) dmMesh, &name);CHECK_PETSC_ERROR(err);
+        msg << "Inconsistency in topology found for mesh '" << name << "' during output.\n"
+            << "Number of vertices (" << coneSize << ") in cell '" << c
+            << "' does not expected number of vertices (" << numCorners << ").";
+        throw std::runtime_error(msg.str());
+      } // if
+      for(PetscInt corner = 0; corner < numCorners; ++corner) {
+        tmpVertices[k++] = vertices[cone[corner]-vStart];
+      }
+    }
+    err = ISRestoreIndices(globalCellNumbers, &cells);CHECK_PETSC_ERROR(err);
+    err = ISRestoreIndices(globalVertexNumbers, &vertices);CHECK_PETSC_ERROR(err);
+#endif
+
     PetscVec elemVec;
     err = VecCreateMPIWithArray(((PetscObject) dmMesh)->comm, numCorners, conesSize, PETSC_DETERMINE,
-				tmpVertices, &elemVec); CHECK_PETSC_ERROR(err);
-    err = PetscObjectSetName((PetscObject) elemVec,
-			     "cells");CHECK_PETSC_ERROR(err);
+                                tmpVertices, &elemVec); CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject) elemVec, "cells");CHECK_PETSC_ERROR(err);
 
     const std::string& filenameCells = _datasetFilename("cells");
     err = PetscViewerBinaryOpen(((PetscObject) dmMesh)->comm, filenameCells.c_str(),
-				FILE_MODE_WRITE,
-				&binaryViewer);
-    CHECK_PETSC_ERROR(err);
-    err = PetscViewerBinarySetSkipHeader(binaryViewer, PETSC_TRUE);
-    CHECK_PETSC_ERROR(err);
-    err = VecView(elemVec, binaryViewer); CHECK_PETSC_ERROR(err);
-    err = VecDestroy(&elemVec); CHECK_PETSC_ERROR(err);
-    err = PetscFree(tmpVertices); CHECK_PETSC_ERROR(err);
-    err = PetscViewerDestroy(&binaryViewer); CHECK_PETSC_ERROR(err);
+                                FILE_MODE_WRITE, &binaryViewer);CHECK_PETSC_ERROR(err);
+    err = PetscViewerBinarySetSkipHeader(binaryViewer, PETSC_TRUE);CHECK_PETSC_ERROR(err);
+    err = VecView(elemVec, binaryViewer);CHECK_PETSC_ERROR(err);
+    err = VecDestroy(&elemVec);CHECK_PETSC_ERROR(err);
+    err = PetscFree(tmpVertices);CHECK_PETSC_ERROR(err);
+    err = PetscViewerDestroy(&binaryViewer);CHECK_PETSC_ERROR(err);
 
     // Create external dataset for cells
     if (!commRank) {
       const hsize_t ndims = 2;
       hsize_t dims[ndims];
-      dims[0] = cNumbering->getGlobalSize();
+      dims[0] = numCells;
       dims[1] = numCorners;
-      _h5->createDatasetRawExternal("/topology", "cells", filenameCells.c_str(),
-				    dims, ndims, scalartype);
+      _h5->createDatasetRawExternal("/topology", "cells", filenameCells.c_str(), dims, ndims, scalartype);
       const int cellDim = mesh.dimension();
-      _h5->writeAttribute("/topology/cells", "cell_dim", (void*)&cellDim,
-			  H5T_NATIVE_INT);
+      _h5->writeAttribute("/topology/cells", "cell_dim", (void*)&cellDim, H5T_NATIVE_INT);
     } // if
     
   } catch (const std::exception& err) {

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -223,11 +223,10 @@
     if (complexMesh) {
       /* DMComplex */
       PetscContainer c;
-      PetscSection   s;
-      Vec            v;
+      PetscSection   s = field.petscSection();
+      Vec            v = field.localVector();
+      assert(s);assert(v);
 
-      field.dmSection(&s, &v);
-
       /* Will change to just VecView() once I setup the vectors correctly (use VecSetOperation() to change the view method) */
       PetscViewerVTKFieldType ft = field.vectorFieldType() != topology::FieldBase::VECTOR ? PETSC_VTK_POINT_FIELD : PETSC_VTK_POINT_VECTOR_FIELD;
       PetscErrorCode err = PetscViewerVTKAddField(_viewer, (PetscObject) complexMesh, DMComplexVTKWriteAll, ft, (PetscObject) v);CHECK_PETSC_ERROR(err);
@@ -236,8 +235,6 @@
       err = PetscContainerSetPointer(c, s);CHECK_PETSC_ERROR(err);
       err = PetscObjectCompose((PetscObject) v, "section", (PetscObject) c);CHECK_PETSC_ERROR(err);
       err = PetscContainerDestroy(&c);CHECK_PETSC_ERROR(err);
-      //err = PetscSectionDestroy(&s);CHECK_PETSC_ERROR(err);
-      err = VecDestroy(&v);CHECK_PETSC_ERROR(err);
 
       _wroteVertexHeader = true;
     } else {
@@ -313,11 +310,10 @@
     if (complexMesh) {
       /* DMComplex */
       PetscContainer c;
-      PetscSection s;
-      Vec v;
+      PetscSection   s = field.petscSection();
+      Vec            v = field.localVector();
+      assert(s);assert(v);
 
-      field.dmSection(&s, &v);
-
       /* Will change to just VecView() once I setup the vectors correctly (use VecSetOperation() to change the view) */
       PetscViewerVTKFieldType ft = field.vectorFieldType() != topology::FieldBase::VECTOR ? PETSC_VTK_CELL_FIELD : PETSC_VTK_CELL_VECTOR_FIELD;
       PetscErrorCode err = PetscViewerVTKAddField(_viewer, (PetscObject) complexMesh, DMComplexVTKWriteAll, ft, (PetscObject) v); CHECK_PETSC_ERROR(err);
@@ -326,8 +322,6 @@
       err = PetscContainerSetPointer(c, s);CHECK_PETSC_ERROR(err);
       err = PetscObjectCompose((PetscObject) v, "section", (PetscObject) c);CHECK_PETSC_ERROR(err);
       err = PetscContainerDestroy(&c);CHECK_PETSC_ERROR(err);
-      //err = PetscSectionDestroy(&s);CHECK_PETSC_ERROR(err);
-      err = VecDestroy(&v);CHECK_PETSC_ERROR(err);
 
       _wroteCellHeader = true;
     } else {

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -198,13 +198,13 @@
   PetscInt     coordSize;
 
   err = DMComplexGetCoordinateSection(complexMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMComplexGetCoordinateVec(complexMesh, &coordVec);CHECK_PETSC_ERROR(err);
   err = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHECK_PETSC_ERROR(err);
   for(PetscInt v = numCells; v < numCells+numVertices; ++v) {
     err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
   }
   err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
   err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+  err = VecCreate(comm, &coordVec);CHECK_PETSC_ERROR(err);
   err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
   err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
   err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
@@ -217,6 +217,7 @@
     }
   }
   err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+  err = DMSetCoordinatesLocal(complexMesh, coordVec);CHECK_PETSC_ERROR(err);
   logger.stagePop(); // Coordinates
 
   sieveMesh->getFactory()->clear();

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIO.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshIO.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -330,8 +330,14 @@
   } // if/else
   sieveMesh->allocate(groupField);
 
-  for(PetscInt p = 0; p < numPoints; ++p) {
-    err = DMComplexSetLabelValue(complexMesh, name.c_str(), points[p], 1);CHECK_PETSC_ERROR(err);
+  if (CELL == type) {
+    for(PetscInt p = 0; p < numPoints; ++p) {
+      err = DMComplexSetLabelValue(complexMesh, name.c_str(), points[p], 1);CHECK_PETSC_ERROR(err);
+    }
+  } else if (VERTEX == type) {
+    for(PetscInt p = 0; p < numPoints; ++p) {
+      err = DMComplexSetLabelValue(complexMesh, name.c_str(), numCells+points[p], 1);CHECK_PETSC_ERROR(err);
+    }
   }
   logger.stagePop();
 } // _setGroup

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/VertexFilterVecNorm.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/VertexFilterVecNorm.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/VertexFilterVecNorm.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -76,18 +76,19 @@
   typedef typename field_type::Mesh::SieveMesh SieveMesh;
   typedef typename SieveMesh::label_sequence label_sequence;
 
-  const ALE::Obj<SieveMesh>& sieveMesh = fieldIn.mesh().sieveMesh();
-  assert(!sieveMesh.isNull());
+  DM dmMesh = fieldIn.mesh().dmMesh();
+  PetscInt       vStart, vEnd;
+  PetscErrorCode err;
 
-  const ALE::Obj<label_sequence>& vertices = sieveMesh->depthStratum(0);
-  assert(!vertices.isNull());
-  const typename label_sequence::iterator verticesEnd = vertices->end();
+  assert(dmMesh);
+  err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
 
-  const ALE::Obj<RealSection>& sectionIn = fieldIn.section();
-  assert(!sectionIn.isNull());
-  const int fiberDimIn = (vertices->size() > 0) ? 
-    sectionIn->getFiberDimension(*vertices->begin()) : 0;
-  const int fiberDimNorm = 1;
+  PetscSection sectionIn = fieldIn.petscSection();
+  Vec          vecIn     = fieldIn.localVector();
+  PetscScalar *a, *an;
+  PetscInt     fiberDimIn, fiberDimNorm = 1;
+  assert(sectionIn);assert(vecIn);
+  err = PetscSectionGetDof(sectionIn, vStart, &fiberDimIn);CHECK_PETSC_ERROR(err);
 
   // Allocation field if necessary
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
@@ -121,25 +122,28 @@
   } // if
   logger.stagePop();
 
-  const ALE::Obj<RealSection>& sectionNorm = 
-    _fieldVecNorm->section();
-  assert(!sectionNorm.isNull());
+  PetscSection sectionNorm = _fieldVecNorm->petscSection();
+  Vec          vecNorm     = _fieldVecNorm->localVector();
+  assert(sectionNorm);assert(vecNorm);
 
-  PylithScalar norm = 0.0;
   // Loop over vertices
-  for (typename label_sequence::iterator v_iter=vertices->begin();
-       v_iter != verticesEnd;
-       ++v_iter) {
-    const PylithScalar* values = sectionIn->restrictPoint(*v_iter);
-    
-    norm = 0.0;
-    for (int i=0; i < fiberDimIn; ++i)
-      norm += values[i]*values[i];
+  err = VecGetArray(vecIn, &a);CHECK_PETSC_ERROR(err);
+  err = VecGetArray(vecNorm, &an);CHECK_PETSC_ERROR(err);
+  for(PetscInt v = vStart; v < vEnd; ++v) {
+    PetscInt dof, off, offn;
+    PylithScalar norm = 0.0;
+
+    err = PetscSectionGetDof(sectionIn, v, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(sectionIn, v, &off);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(sectionNorm, v, &offn);CHECK_PETSC_ERROR(err);
+    for(PetscInt d = 0; d < dof; ++d)
+      norm += a[off+d]*a[off+d];
     norm = sqrt(norm);
-
-    sectionNorm->updatePoint(*v_iter, &norm);
+    an[offn] = norm;
   } // for
-  PetscLogFlops(vertices->size() * (1 + 2*fiberDimIn) );
+  err = VecRestoreArray(vecIn, &a);CHECK_PETSC_ERROR(err);
+  err = VecRestoreArray(vecNorm, &an);CHECK_PETSC_ERROR(err);
+  PetscLogFlops((vEnd-vStart) * (1 + 2*fiberDimIn));
 
   return *_fieldVecNorm;
 } // filter

Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -295,7 +295,7 @@
 
     err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
     err = DMComplexGetCoordinateSection(dmMesh, &coordinateSection);CHECK_PETSC_ERROR(err);
-    err = DMComplexGetCoordinateVec(dmMesh, &coordinateVec);CHECK_PETSC_ERROR(err);
+    err = DMGetCoordinatesLocal(dmMesh, &coordinateVec);CHECK_PETSC_ERROR(err);
     assert(coordinateSection);assert(coordinateVec);
     if (dim > 1) {
       const int m = (dim * (dim + 1)) / 2;

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -64,17 +64,16 @@
   assert(!section.isNull());
   _metadata["default"] = metadata;
   if (mesh.dmMesh()) {
-    Vec          lv;
     PetscSection s;
     PetscErrorCode err = DMComplexClone(mesh.dmMesh(), &_dm);CHECK_PETSC_ERROR(err);
-    this->dmSection(&s, &lv);
-    err = VecDestroy(&lv);CHECK_PETSC_ERROR(err);
+    this->dmSection(&s, &_localVec);
     err = DMSetDefaultSection(_dm, s);CHECK_PETSC_ERROR(err);
+    err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject) _localVec,  _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
   } else {
     _dm = PETSC_NULL;
   }
-  _globalVec = PETSC_NULL;
-  _localVec  = PETSC_NULL;
 } // constructor
 
 // ----------------------------------------------------------------------
@@ -143,6 +142,12 @@
   if (!_section.isNull()) {
     _section->setName(value);
   } // if
+  if (_localVec) {
+    PetscErrorCode err = PetscObjectSetName((PetscObject) _localVec, value);CHECK_PETSC_ERROR(err);
+  }
+  if (_globalVec) {
+    PetscErrorCode err = PetscObjectSetName((PetscObject) _globalVec, value);CHECK_PETSC_ERROR(err);
+  }
 
   const typename scatter_map_type::const_iterator scattersEnd = _scatters.end();
   for (typename scatter_map_type::const_iterator s_iter=_scatters.begin();
@@ -325,6 +330,33 @@
 					       const int fiberDim,
 					       const int stratum)
 { // newSection
+  // Changing this because cells/vertices are numbered differently in the new scheme
+  if (_dm) {
+    PetscSection   s;
+    PetscInt       pStart, pEnd;
+    PetscErrorCode err;
+
+    switch(domain) {
+    case VERTICES_FIELD:
+      err = DMComplexGetDepthStratum(_dm, stratum, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+      break;
+    case CELLS_FIELD:
+      err = DMComplexGetHeightStratum(_dm, stratum, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+      break;
+    case POINTS_FIELD:
+      err = DMComplexGetChart(_dm, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+      break;
+    default:
+      std::cerr << "Unknown value for DomainEnum: " << domain << std::endl;
+      throw std::logic_error("Bad domain enum in Field.");
+    }
+    err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+    err = PetscSectionSetChart(s, pStart, pEnd);CHECK_PETSC_ERROR(err);
+
+    for(PetscInt p = pStart; p < pEnd; ++p) {
+      err = PetscSectionSetDof(s, p, fiberDim);CHECK_PETSC_ERROR(err);
+    }
+  } else {
   const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
   assert(!sieveMesh.isNull());
 
@@ -340,6 +372,7 @@
   } // else
 
   newSection(points, fiberDim);
+  }
 } // newSection
 
 // ----------------------------------------------------------------------
@@ -391,6 +424,18 @@
         err = PetscSectionSetDof(s, p, fiberDim);CHECK_PETSC_ERROR(err);
       }
     }
+  } else if (src._dm) {
+    PetscSection   srcs, s;
+    PetscInt       pStart, pEnd;
+    PetscErrorCode err;
+
+    err = DMGetDefaultSection(src._dm, &srcs);CHECK_PETSC_ERROR(err);
+    err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetChart(srcs, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+    err = PetscSectionSetChart(s, pStart, pEnd);CHECK_PETSC_ERROR(err);
+    for(PetscInt p = pStart; p < pEnd; ++p) {
+      err = PetscSectionSetDof(s, p, fiberDim);CHECK_PETSC_ERROR(err);
+    }
   } // if
 
   logger.stagePop();
@@ -442,16 +487,19 @@
       _section->setBC(srcSection->getBC());
       _section->copySpaces(srcSection); // :BUG: NEED TO REBUILD SPACES 
     } // if/else
-    PetscSection   section = src.petscSection();
-    PetscSection   newSection;
-    PetscErrorCode err;
+  } // if
+  PetscSection   section = src.petscSection();
+  PetscSection   newSection;
+  PetscErrorCode err;
 
-    if (_dm) {
-      err = PetscSectionClone(section, &newSection);CHECK_PETSC_ERROR(err);
-      err = DMSetDefaultSection(_dm, newSection);CHECK_PETSC_ERROR(err);
-      err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
-      err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
-    }
+  if (_dm) {
+    err = PetscSectionClone(section, &newSection);CHECK_PETSC_ERROR(err);
+    err = DMSetDefaultSection(_dm, newSection);CHECK_PETSC_ERROR(err);
+    err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
+    err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject) _localVec,  _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+  }
     
     // Reuse scatters in clone
     if (src._scatters.size() > 0) {
@@ -512,7 +560,6 @@
 	_scatters[s_iter->first] = sinfo;
       } // for
     } // if
-  } // if
   logger.stagePop();
 } // cloneSection
 
@@ -688,6 +735,8 @@
     err = PetscSectionSetUp(s);CHECK_PETSC_ERROR(err);
     err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
     err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject) _localVec,  _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
   }
 
   logger.stagePop();
@@ -1044,10 +1093,8 @@
 								const char* context)
 { // createScatter
   assert(context);
-  assert(!_section.isNull());
-  assert(!mesh.sieveMesh().isNull());
+  PetscErrorCode err = 0;
 
-  PetscErrorCode err = 0;
   const bool createScatterOk = true;
   ScatterInfo& sinfo = _getScatter(context, createScatterOk);
   if (sinfo.scatter) {
@@ -1059,45 +1106,42 @@
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("GlobalOrder");
 
-  // Get global order (create if necessary).
-  const std::string& orderLabel = _section->getName();
-  const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
-    mesh.sieveMesh();
-  assert(!sieveMesh.isNull());
-  const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order = 
-    sieveMesh->getFactory()->getGlobalOrder(sieveMesh, orderLabel,
-					    _section);
-  assert(!order.isNull());
+  if (!_section.isNull()) {
+    assert(!mesh.sieveMesh().isNull());
+    // Get global order (create if necessary).
+    const std::string& orderLabel = _section->getName();
+    const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
+    assert(!sieveMesh.isNull());
+    const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order = 
+      sieveMesh->getFactory()->getGlobalOrder(sieveMesh, orderLabel, _section);
+    assert(!order.isNull());
 
-  // Create scatter
-  err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, false, 
-				  &sinfo.scatter);
-  CHECK_PETSC_ERROR(err);
+    // Create scatter
+    err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, false, &sinfo.scatter);
+    CHECK_PETSC_ERROR(err);
   
-  // Create scatterVec
-  const int blockSize = 1;
-  if (_section->sizeWithBC() > 0) {
-    err = VecCreateSeqWithArray(PETSC_COMM_SELF,
-				blockSize, _section->getStorageSize(),
-				_section->restrictSpace(),
-				&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
-  } else {
-    err = VecCreateSeqWithArray(PETSC_COMM_SELF, 
-				blockSize, 0, PETSC_NULL,
-				&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
-  } // else
+    // Create scatterVec
+    const int blockSize = 1;
+    if (_section->sizeWithBC() > 0) {
+      err = VecCreateSeqWithArray(PETSC_COMM_SELF,
+                                  blockSize, _section->getStorageSize(),
+                                  _section->restrictSpace(),
+                                  &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } else {
+      err = VecCreateSeqWithArray(PETSC_COMM_SELF, 
+                                  blockSize, 0, PETSC_NULL,
+                                  &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } // else
 
-  // Create vector
 #if 0
-  err = VecCreate(_mesh.comm(), &sinfo.vector);
-  CHECK_PETSC_ERROR(err);
-  err = PetscObjectSetName((PetscObject)sinfo.vector,
-			   _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(sinfo.vector,
-		    order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
-  err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(sinfo.vector);CHECK_PETSC_ERROR(err);
+    // Create vector
+    err = VecCreate(_mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+    err = VecSetSizes(sinfo.vector, order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
+    err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
+    err = VecSetFromOptions(sinfo.vector);CHECK_PETSC_ERROR(err);
 #endif
+  }
   PetscInt localSize, globalSize;
 
   err = PetscObjectReference((PetscObject) _dm);CHECK_PETSC_ERROR(err);
@@ -1105,7 +1149,7 @@
   err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
   err = VecGetSize(_localVec,  &localSize);CHECK_PETSC_ERROR(err);
   err = VecGetSize(_globalVec, &globalSize);CHECK_PETSC_ERROR(err);
-  assert(order->getLocalSize()  == localSize);
+  //assert(order->getLocalSize()  == localSize);
   //assert(order->getGlobalSize() == globalSize);
   sinfo.vector = _globalVec;
   sinfo.dm     = _dm;
@@ -1129,10 +1173,8 @@
 { // createScatter
   assert(!numbering.isNull());
   assert(context);
-  assert(!_section.isNull());
-  assert(!mesh.sieveMesh().isNull());
+  PetscErrorCode err = 0;
 
-  PetscErrorCode err = 0;
   const bool createScatterOk = true;
   ScatterInfo& sinfo = _getScatter(context, createScatterOk);
   
@@ -1146,49 +1188,47 @@
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("GlobalOrder");
 
-  // Get global order (create if necessary).
-  const std::string& orderLabel = 
-    (strlen(context) > 0) ?
-    _section->getName() + std::string("_") + std::string(context) :
-    _section->getName();
-  const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
-    mesh.sieveMesh();
-  assert(!sieveMesh.isNull());
-  const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order = 
-    sieveMesh->getFactory()->getGlobalOrder(sieveMesh, orderLabel,
-                                            numbering->getChart().begin(),
-                                            numbering->getChart().end(),
-                                            _section);
-  assert(!order.isNull());
+  if (!_section.isNull()); {
+    assert(!mesh.sieveMesh().isNull());
+    // Get global order (create if necessary).
+    const std::string& orderLabel = 
+      (strlen(context) > 0) ?
+      _section->getName() + std::string("_") + std::string(context) :
+      _section->getName();
+    const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
+    assert(!sieveMesh.isNull());
+    const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order = 
+      sieveMesh->getFactory()->getGlobalOrder(sieveMesh, orderLabel,
+                                              numbering->getChart().begin(),
+                                              numbering->getChart().end(),
+                                              _section);
+    assert(!order.isNull());
 
-  // Create scatter
-  err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, false, 
-				  &sinfo.scatter); 
-  CHECK_PETSC_ERROR(err);
+    // Create scatter
+    err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, false, &sinfo.scatter);CHECK_PETSC_ERROR(err);
 
-  // Create scatterVec
-  const int blockSize = 1;
-  if (_section->sizeWithBC() > 0) {
-    err = VecCreateSeqWithArray(PETSC_COMM_SELF,
-				blockSize, _section->getStorageSize(),
-				_section->restrictSpace(),
-				&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
-  } else {
-    err = VecCreateSeqWithArray(PETSC_COMM_SELF, 
-				blockSize, 0, PETSC_NULL,
-				&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
-  } // else
+    // Create scatterVec
+    const int blockSize = 1;
+    if (_section->sizeWithBC() > 0) {
+      err = VecCreateSeqWithArray(PETSC_COMM_SELF,
+                                  blockSize, _section->getStorageSize(),
+                                  _section->restrictSpace(),
+                                  &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } else {
+      err = VecCreateSeqWithArray(PETSC_COMM_SELF, 
+                                  blockSize, 0, PETSC_NULL,
+                                  &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } // else
 
-  // Create vector
 #if 0
-  err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
-  err = PetscObjectSetName((PetscObject)sinfo.vector,
-			   _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(sinfo.vector,
-		    order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
-  err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);  
+    // Create vector
+    err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+    err = VecSetSizes(sinfo.vector, order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
+    err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
+    err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);  
 #endif
+  }
   PetscInt localSize, globalSize;
 
   err = PetscObjectReference((PetscObject) _dm);CHECK_PETSC_ERROR(err);
@@ -1196,8 +1236,8 @@
   err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
   err = VecGetSize(_localVec,  &localSize);CHECK_PETSC_ERROR(err);
   err = VecGetSize(_globalVec, &globalSize);CHECK_PETSC_ERROR(err);
-  assert(order->getLocalSize()  == localSize);
-  assert(order->getGlobalSize() == globalSize);
+  //assert(order->getLocalSize()  == localSize);
+  //assert(order->getGlobalSize() == globalSize);
   sinfo.vector = _globalVec;
   sinfo.dm     = _dm;
 
@@ -1228,11 +1268,8 @@
 	const char* context)
 { // createScatterWithBC
   assert(context);
-  assert(!_section.isNull());
-  assert(!mesh.sieveMesh().isNull());
+  PetscErrorCode err = 0;
 
-
-  PetscErrorCode err = 0;
   const bool createScatterOk = true;
   ScatterInfo& sinfo = _getScatter(context, createScatterOk);
   if (sinfo.scatter) {
@@ -1244,42 +1281,42 @@
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("GlobalOrder");
 
-  // Get global order (create if necessary).
-  const std::string& orderLabel = _section->getName();
-  const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
-    mesh.sieveMesh();
-  assert(!sieveMesh.isNull());
-  const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order = 
-    sieveMesh->getFactory()->getGlobalOrderWithBC(sieveMesh, orderLabel,
-						  _section);
-  assert(!order.isNull());
+  if (!_section.isNull()) {
+    assert(!mesh.sieveMesh().isNull());
 
-  // Create scatter
-  err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, true, 
-				  &sinfo.scatter); 
-  CHECK_PETSC_ERROR(err);
+    // Get global order (create if necessary).
+    const std::string& orderLabel = _section->getName();
+    const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
+    assert(!sieveMesh.isNull());
+    const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order = 
+      sieveMesh->getFactory()->getGlobalOrderWithBC(sieveMesh, orderLabel, _section);
+    assert(!order.isNull());
+
+    // Create scatter
+    err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, true, &sinfo.scatter); 
+    CHECK_PETSC_ERROR(err);
   
-  // Create scatterVec
-  const int blockSize = _getFiberDim();
-  if (_section->sizeWithBC() > 0) {
-    err = VecCreateSeqWithArray(PETSC_COMM_SELF,
-				blockSize, _section->getStorageSize(),
-				_section->restrictSpace(),
-				&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
-  } else {
-    err = VecCreateSeqWithArray(PETSC_COMM_SELF, 
-				blockSize, 0, PETSC_NULL,
-				&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
-  } // else
-  
-  // Create vector
+    // Create scatterVec
+    const int blockSize = _getFiberDim();
+    if (_section->sizeWithBC() > 0) {
+      err = VecCreateSeqWithArray(PETSC_COMM_SELF,
+                                  blockSize, _section->getStorageSize(),
+                                  _section->restrictSpace(),
+                                  &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } else {
+      err = VecCreateSeqWithArray(PETSC_COMM_SELF, 
+                                  blockSize, 0, PETSC_NULL,
+                                  &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } // else
 #if 0
-  err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
-  err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(sinfo.vector, order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
-  err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(sinfo.vector);CHECK_PETSC_ERROR(err);
+    // Create vector
+    err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+    err = VecSetSizes(sinfo.vector, order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
+    err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
+    err = VecSetFromOptions(sinfo.vector);CHECK_PETSC_ERROR(err);
 #endif
+  }
 
   PetscSection section, newSection, gsection;
   PetscSF      sf;
@@ -1297,8 +1334,8 @@
 
   err = PetscSectionGetStorageSize(section, &localSize);CHECK_PETSC_ERROR(err);
   err = VecGetSize(sinfo.vector, &globalSize);CHECK_PETSC_ERROR(err);
-  assert(order->getLocalSize()  == localSize);
-  assert(order->getGlobalSize() == globalSize);
+  //assert(order->getLocalSize()  == localSize);
+  //assert(order->getGlobalSize() == globalSize);
 
   logger.stagePop();
 } // createScatterWithBC
@@ -1319,9 +1356,8 @@
 { // createScatterWithBC
   assert(!numbering.isNull());
   assert(context);
-  assert(!_section.isNull());
+  PetscErrorCode err = 0;
 
-  PetscErrorCode err = 0;
   const bool createScatterOk = true;
   ScatterInfo& sinfo = _getScatter(context, createScatterOk);
   
@@ -1335,58 +1371,64 @@
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("GlobalOrder");
 
-  // Get global order (create if necessary).
-  const std::string& orderLabel = 
-    (strlen(context) > 0) ?
-    _section->getName() + std::string("_") + std::string(context) :
-    _section->getName();
-  const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
-    mesh.sieveMesh();
-  assert(!sieveMesh.isNull());
-  const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order = 
-    sieveMesh->getFactory()->getGlobalOrderWithBC(sieveMesh, orderLabel,
-                                                  numbering->getChart().begin(),
-                                                  numbering->getChart().end(),
-                                                  _section);
-  assert(!order.isNull());
-  //order->view("GLOBAL ORDER"); // DEBUG
+  if (!_section.isNull()) {
+    // Get global order (create if necessary).
+    const std::string& orderLabel = 
+      (strlen(context) > 0) ?
+      _section->getName() + std::string("_") + std::string(context) :
+      _section->getName();
+    const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
+    assert(!sieveMesh.isNull());
+    const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order = 
+      sieveMesh->getFactory()->getGlobalOrderWithBC(sieveMesh, orderLabel,
+                                                    numbering->getChart().begin(),
+                                                    numbering->getChart().end(),
+                                                    _section);
+    assert(!order.isNull());
+    //order->view("GLOBAL ORDER"); // DEBUG
 
-  // Create scatter
-  err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, true, 
-				  &sinfo.scatter); 
-  CHECK_PETSC_ERROR(err);
+    // Create scatter
+    err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, true, &sinfo.scatter); 
+    CHECK_PETSC_ERROR(err);
 
-  // Create scatterVec
-  const int blockSize = _getFiberDim();
-  if (_section->sizeWithBC() > 0) {
-    err = VecCreateSeqWithArray(PETSC_COMM_SELF,
-				blockSize, _section->getStorageSize(),
-				_section->restrictSpace(),
-				&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
-  } else {
-    err = VecCreateSeqWithArray(PETSC_COMM_SELF, 
-				blockSize, 0, PETSC_NULL,
-				&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
-  } // else
+    // Create scatterVec
+    const int blockSize = _getFiberDim();
+    if (_section->sizeWithBC() > 0) {
+      err = VecCreateSeqWithArray(PETSC_COMM_SELF,
+                                  blockSize, _section->getStorageSize(),
+                                  _section->restrictSpace(),
+                                  &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } else {
+      err = VecCreateSeqWithArray(PETSC_COMM_SELF, 
+                                  blockSize, 0, PETSC_NULL,
+                                  &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } // else
 
-  // Create vector
 #if 0
-  err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
-  err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
-  err = VecSetSizes(sinfo.vector,order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
-  err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
-  err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);  
+    // Create vector
+    err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
+    err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+    err = VecSetSizes(sinfo.vector,order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
+    err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
+    err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);  
 #endif
+  }
 
   PetscSection section, newSection, gsection;
   PetscSF      sf;
+  PetscInt     cEnd, cMax, vEnd, vMax;
+  err = DMComplexGetHeightStratum(_dm, 0, PETSC_NULL, &cEnd);CHECK_PETSC_ERROR(err);
+  err = DMComplexGetDepthStratum(_dm, 0, PETSC_NULL, &vEnd);CHECK_PETSC_ERROR(err);
+  err = DMComplexGetVTKBounds(_dm, &cMax, &vMax);CHECK_PETSC_ERROR(err);
+  PetscInt     excludeRanges[4] = {cMax, cEnd, vMax, vEnd};
+  PetscInt     numExcludes      = (cMax >= 0 ? 1 : 0) + (vMax >= 0 ? 1 : 0);
 
   err = DMComplexClone(_dm, &sinfo.dm);CHECK_PETSC_ERROR(err);
   err = DMGetDefaultSection(_dm, &section);CHECK_PETSC_ERROR(err);
   err = PetscSectionClone(section, &newSection);CHECK_PETSC_ERROR(err);
   err = DMSetDefaultSection(sinfo.dm, newSection);CHECK_PETSC_ERROR(err);
   err = DMGetPointSF(sinfo.dm, &sf);CHECK_PETSC_ERROR(err);
-  err = PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, &gsection);CHECK_PETSC_ERROR(err);
+  err = PetscSectionCreateGlobalSectionCensored(section, sf, PETSC_TRUE, numExcludes, excludeRanges, &gsection);CHECK_PETSC_ERROR(err);
   err = DMSetDefaultGlobalSection(sinfo.dm, gsection);CHECK_PETSC_ERROR(err);
   err = DMCreateGlobalVector(sinfo.dm, &sinfo.vector);CHECK_PETSC_ERROR(err);
   err = PetscObjectSetName((PetscObject) sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
@@ -1394,8 +1436,8 @@
 
   err = PetscSectionGetStorageSize(section, &localSize);CHECK_PETSC_ERROR(err);
   err = VecGetSize(sinfo.vector, &globalSize);CHECK_PETSC_ERROR(err);
-  assert(order->getLocalSize()  == localSize);
-  assert(order->getGlobalSize() == globalSize);
+  /* assert(order->getLocalSize()  == localSize); This does not work because the local vector includes the lagrange cell variables */
+  /* assert(order->getGlobalSize() == globalSize); */
 
 #if 0
   std::cout << "["<<sieveMesh->commRank()<<"] CONTEXT: " << context 
@@ -1462,17 +1504,16 @@
 { // scatterSectionToVector
   assert(vector);
   assert(context);
-  assert(!_section.isNull());
-
-  PetscErrorCode err = 0;
   const ScatterInfo& sinfo = _getScatter(context);
+  PetscErrorCode     err   = 0;
 #if 0
-  err = VecScatterBegin(sinfo.scatter, sinfo.scatterVec, vector,
-			INSERT_VALUES, SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
-  err = VecScatterEnd(sinfo.scatter, sinfo.scatterVec, vector,
-		      INSERT_VALUES, SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
+  if (!_section.isNull()) {
+    err = VecScatterBegin(sinfo.scatter, sinfo.scatterVec, vector,
+                          INSERT_VALUES, SCATTER_FORWARD);CHECK_PETSC_ERROR(err);
+    err = VecScatterEnd(sinfo.scatter, sinfo.scatterVec, vector,
+                        INSERT_VALUES, SCATTER_FORWARD);CHECK_PETSC_ERROR(err);
+  }
 #endif
-
   if (sinfo.dm) {
     err = DMLocalToGlobalBegin(sinfo.dm, _localVec, INSERT_VALUES, vector);CHECK_PETSC_ERROR(err);
     err = DMLocalToGlobalEnd(sinfo.dm, _localVec, INSERT_VALUES, vector);CHECK_PETSC_ERROR(err);
@@ -1502,17 +1543,17 @@
 { // scatterVectorToSection
   assert(vector);
   assert(context);
-  assert(!_section.isNull());
+  const ScatterInfo& sinfo = _getScatter(context);
+  PetscErrorCode     err   = 0;
 
-  PetscErrorCode err = 0;
-  const ScatterInfo& sinfo = _getScatter(context);
 #if 0
-  err = VecScatterBegin(sinfo.scatter, vector, sinfo.scatterVec,
-			INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
-  err = VecScatterEnd(sinfo.scatter, vector, sinfo.scatterVec,
-		      INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
+  if (!_section.isNull()) {
+    err = VecScatterBegin(sinfo.scatter, vector, sinfo.scatterVec,
+                          INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
+    err = VecScatterEnd(sinfo.scatter, vector, sinfo.scatterVec,
+                        INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
+  }
 #endif
-
   if (sinfo.dm) {
     err = DMGlobalToLocalBegin(sinfo.dm, vector, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
     err = DMGlobalToLocalEnd(sinfo.dm, vector, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
@@ -1685,6 +1726,9 @@
   case CELLS_FIELD:
     err = DMComplexGetHeightStratum(_dm, 0, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
     break;
+  case POINTS_FIELD:
+    err = DMComplexGetChart(_dm, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+    break;
   default:
     std::cerr << "Unknown value for DomainEnum: " << domain << std::endl;
     throw std::logic_error("Bad domain enum in Field.");

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/FieldBase.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/FieldBase.hh	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/FieldBase.hh	2012-10-04 12:17:00 UTC (rev 20803)
@@ -55,6 +55,7 @@
   enum DomainEnum {
     VERTICES_FIELD=0, ///< FieldBase over vertices.
     CELLS_FIELD=1, ///< FieldBase over cells.
+    POINTS_FIELD=2, ///< FieldBase over all points.
   }; // DomainEnum
 
 // PUBLIC STRUCTS ///////////////////////////////////////////////////////

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -20,6 +20,8 @@
 #error "Mesh.icc must be included only from Mesh.hh"
 #else
 
+#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
+
 // Get Sieve mesh.
 inline
 const ALE::Obj<pylith::topology::Mesh::SieveMesh>&
@@ -45,7 +47,10 @@
 inline
 void
 pylith::topology::Mesh::setDMMesh(DM dm) {
+  PetscErrorCode err;
+  err = DMDestroy(&_newMesh);CHECK_PETSC_ERROR(err);
   _newMesh = dm;
+  err = PetscObjectSetName((PetscObject) _newMesh, "domain");CHECK_PETSC_ERROR(err);
 }
 
 // Get point type sizes.

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -86,6 +86,7 @@
   } // if
   PetscBool      hasLabel;
   PetscErrorCode err;
+
   err = DMComplexHasLabel(dmMesh, label, &hasLabel);CHECK_PETSC_ERROR(err);
   if (!hasLabel) {
     std::ostringstream msg;
@@ -110,8 +111,8 @@
     _mesh->setRealSection("coordinates_dimensioned", 
 			  meshSieveMesh->getRealSection("coordinates_dimensioned"));
 
-  /* TODO: Implement subMesh() for DMComplex */
-  _newMesh = PETSC_NULL;
+  /* TODO: Add creation of pointSF for submesh */
+  err = DMComplexCreateSubmesh(dmMesh, label, &_newMesh);CHECK_PETSC_ERROR(err);
 
   // Create the parallel overlap
   const ALE::Obj<SieveMesh::sieve_type>& sieve = _mesh->getSieve();
@@ -155,6 +156,7 @@
   // Set name
   std::string meshLabel = "subdomain_" + std::string(label);
   _mesh->setName(meshLabel);
+  err = PetscObjectSetName((PetscObject) _newMesh, meshLabel.c_str());CHECK_PETSC_ERROR(err);
 
   int maxConeSizeLocal = sieve->getMaxConeSize();
   int maxConeSize = 0;

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh	2012-10-04 12:17:00 UTC (rev 20803)
@@ -102,6 +102,12 @@
    */
   DM dmMesh(void) const;
 
+  /** Set DMComplex mesh.
+   *
+   * @param DMComplex mesh.
+   */
+  void setDMMesh(DM dm);
+
   /** Get sizes for all point types.
    *
    * @param numNormalCells

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -20,6 +20,8 @@
 #error "SubMesh.icc must be included only from SubMesh.hh"
 #else
 
+#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
+
 // Get Sieve mesh.
 inline
 const ALE::Obj<pylith::topology::Mesh::SieveSubMesh>&
@@ -41,6 +43,16 @@
   return _newMesh;
 }
 
+// Set DMComplex mesh.
+inline
+void
+pylith::topology::SubMesh::setDMMesh(DM dm) {
+  PetscErrorCode err;
+  err = DMDestroy(&_newMesh);CHECK_PETSC_ERROR(err);
+  _newMesh = dm;
+  err = PetscObjectSetName((PetscObject) _newMesh, "domain");CHECK_PETSC_ERROR(err);
+}
+
 // Get point type sizes.
 inline
 void

Modified: short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -155,6 +155,7 @@
   const int fieldsFiberDim = numProperties;
   int index = 0;
   CPPUNIT_ASSERT(0 != friction._fieldsPropsStateVars);
+#if 1
   const ALE::Obj<SubRealUniformSection>& fieldsSection =
     friction._fieldsPropsStateVars->section();
   CPPUNIT_ASSERT(!fieldsSection.isNull());
@@ -173,7 +174,33 @@
         CPPUNIT_ASSERT_DOUBLES_EQUAL(propertiesE[index], fieldsVertex[i],
 				     tolerance);
   } // for
+#else
+  PetscSection fieldsSection = friction._fieldsPropsStateVars->petscSection();
+  Vec          fieldsVec     = friction._fieldsPropsStateVars->localVector();
+  PetscScalar *fieldsArray;
+  PetscErrorCode err;
 
+  CPPUNIT_ASSERT(fieldsSection);CPPUNIT_ASSERT(fieldsVec);
+  err = VecGetArray(fieldsVec, &fieldsArray);CHECK_PETSC_ERROR(err);
+  for(SieveSubMesh::label_sequence::iterator v_iter = verticesBegin;
+      v_iter != verticesEnd;
+      ++v_iter) {
+    PetscInt dof, off;
+
+    err = PetscSectionGetDof(fieldsSection, *v_iter, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(fieldsSection, *v_iter, &off);CHECK_PETSC_ERROR(err);
+    CPPUNIT_ASSERT_EQUAL(fieldsFiberDim, dof);
+    for (int i = 0; i < numProperties; ++i, ++index)
+      if (0 != propertiesE[index])
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, fieldsArray[off+i]/propertiesE[index],
+				     tolerance);
+      else
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(propertiesE[index], fieldsArray[off+i],
+				     tolerance);
+  } // for
+  err = VecRestoreArray(fieldsVec, &fieldsArray);CHECK_PETSC_ERROR(err);
+#endif
+
   // Test vertex array sizes.
   size_t size = data.numPropsVertex + data.numVarsVertex;
   CPPUNIT_ASSERT_EQUAL(size, friction._propsStateVarsVertex.size());
@@ -196,32 +223,34 @@
   const topology::Field<topology::SubMesh>& field = 
     friction.getField("friction_coefficient");
 
-  const ALE::Obj<SieveSubMesh>& sieveMesh = field.mesh().sieveMesh();
-  assert(!sieveMesh.isNull());
-  const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
-    sieveMesh->depthStratum(0);
-  assert(!vertices.isNull());
-  const SieveSubMesh::label_sequence::iterator verticesBegin =
-    vertices->begin();
-  const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+  DM             dmMesh = field.mesh().dmMesh();
+  PetscInt       vStart, vEnd;
+  PetscErrorCode err;
 
+  assert(dmMesh);
+  err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
+  PetscSection fieldSection = field.petscSection();
+  Vec          fieldVec     = field.localVector();
+  PetscScalar *fieldArray;
   const PylithScalar tolerance = 1.0e-06;
+  int index = 0;
 
-  int index = 0;
-  scalar_array fieldVertex(fiberDim);
-  const ALE::Obj<RealSection>& fieldSection = field.section();
-  CPPUNIT_ASSERT(!fieldSection.isNull());
-  for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
-       v_iter != verticesEnd;
-       ++v_iter) {
-    CPPUNIT_ASSERT_EQUAL(fiberDim, fieldSection->getFiberDimension(*v_iter));
-    fieldSection->restrictPoint(*v_iter, &fieldVertex[0], fieldVertex.size());
-    for (int i=0; i < fiberDim; ++i, ++index)
+  CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
+  err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+  for(PetscInt v = vStart; v < vEnd; ++v) {
+    PetscInt dof, off;
+
+    err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
+    CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+    for(int i = 0; i < fiberDim; ++i, ++index)
       if (0 != fieldE[index])
-        CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, fieldVertex[i]/fieldE[index], tolerance);
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, fieldArray[off+i]/fieldE[index], tolerance);
       else
-        CPPUNIT_ASSERT_DOUBLES_EQUAL(fieldE[index], fieldVertex[i], tolerance);
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(fieldE[index], fieldArray[off+i], tolerance);
   } // for
+  err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
 } // testGetField
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestCellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestCellFilterAvg.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestCellFilterAvg.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -100,23 +100,29 @@
   field.vectorFieldType(fieldType);
   field.label(label.c_str());
 
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  const ALE::Obj<SieveMesh::label_sequence>& cells = 
-    sieveMesh->heightStratum(0);
-  CPPUNIT_ASSERT(!cells.isNull());
-  const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+  DM dmMesh = mesh.dmMesh();
+  PetscInt       cStart, cEnd;
+  PetscErrorCode err;
 
-  const ALE::Obj<RealSection>& section = field.section();
-  CPPUNIT_ASSERT(!section.isNull());
-  CPPUNIT_ASSERT_EQUAL(ncells, int(cells->size()));
-  int ipt = 0;
-  for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
-       c_iter != cellsEnd;
-       ++c_iter, ++ipt) {
-    const PylithScalar* values = &fieldValues[ipt*fiberDim];
-    section->updatePoint(*c_iter, values);
+  CPPUNIT_ASSERT(dmMesh);
+  err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
+  PetscSection section = field.petscSection();
+  Vec          vec     = field.localVector();
+  PetscScalar *a;
+  CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+  CPPUNIT_ASSERT_EQUAL(ncells, cEnd-cStart);
+  err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+  for(PetscInt c = cStart; c < cEnd; ++c) {
+    PetscInt dof, off;
+
+    err = PetscSectionGetDof(section, c, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(section, c, &off);CHECK_PETSC_ERROR(err);
+    for(PetscInt d = 0; d < dof; ++d) {
+      a[off+d] = fieldValues[(c-cStart)*dof+d];
+    }
   } // for
+  err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
 
   feassemble::Quadrature<topology::Mesh> quadrature;
   quadrature.initialize(basis, numQuadPts, numBasis,
@@ -129,25 +135,26 @@
   filter.quadrature(&quadrature);
 
   const topology::Field<topology::Mesh>& fieldF = filter.filter(field);
-  const ALE::Obj<RealSection>& sectionF = fieldF.section();
-  CPPUNIT_ASSERT(!sectionF.isNull());
+  PetscSection sectionF = fieldF.petscSection();
+  Vec          vecF     = fieldF.localVector();
+  CPPUNIT_ASSERT(sectionF);CPPUNIT_ASSERT(vecF);
 
   CPPUNIT_ASSERT_EQUAL(fieldTypeE, fieldF.vectorFieldType());
   CPPUNIT_ASSERT_EQUAL(label, std::string(fieldF.label()));
 
-  ipt = 0;
-  for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
-       c_iter != cellsEnd;
-       ++c_iter, ++ipt) {
-    CPPUNIT_ASSERT_EQUAL(fiberDimE, sectionF->getFiberDimension(*c_iter));
-    const PylithScalar* values = sectionF->restrictPoint(*c_iter);
-    CPPUNIT_ASSERT(0 != values);
-    const PylithScalar tolerance = 1.0e-06;
-    for (int i=0; i < fiberDimE; ++i)
-      CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, 
-				   values[i]/fieldValuesE[ipt*fiberDimE+i],
-				   tolerance);
+  const PylithScalar tolerance = 1.0e-06;
+  err = VecGetArray(vecF, &a);CHECK_PETSC_ERROR(err);
+  for(PetscInt c = cStart; c < cEnd; ++c) {
+    PetscInt dof, off;
+
+    err = PetscSectionGetDof(sectionF, c, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(sectionF, c, &off);CHECK_PETSC_ERROR(err);
+    CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
+    for(PetscInt d = 0; d < dof; ++d) {
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, a[off+d]/fieldValuesE[(c-cStart)*dof+d], tolerance);
+    }
   } // for
+  err = VecRestoreArray(vecF, &a);CHECK_PETSC_ERROR(err);
 } // testFilter
 
 

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterBCMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterBCMesh.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterBCMesh.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -110,15 +110,13 @@
   try {
     const int nfields = _data->numVertexFields;
 
-    const ALE::Obj<topology::SubMesh::SieveMesh>& sieveSubMesh = 
-      _submesh->sieveMesh();
-    CPPUNIT_ASSERT(!sieveSubMesh.isNull());
-    const ALE::Obj<topology::SubMesh::SieveMesh::label_sequence>& vertices =
-      sieveSubMesh->depthStratum(0);
-    CPPUNIT_ASSERT(!vertices.isNull());
-    const topology::SubMesh::SieveMesh::label_sequence::iterator verticesEnd =
-      vertices->end();
+    DM dmMesh = _submesh->dmMesh();
+    PetscInt       vStart, vEnd;
+    PetscErrorCode err;
 
+    CPPUNIT_ASSERT(dmMesh);
+    err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
     // Set vertex fields
     for (int i=0; i < nfields; ++i) {
       const char* name = _data->vertexFieldsInfo[i].name;
@@ -129,16 +127,22 @@
       field.allocate();
       field.vectorFieldType(_data->vertexFieldsInfo[i].field_type);
 
-      const ALE::Obj<topology::SubMesh::RealSection>& section = field.section();
-      CPPUNIT_ASSERT(!section.isNull());
-      int ipt = 0;
-      for (topology::SubMesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
-	   v_iter != verticesEnd;
-	   ++v_iter, ++ipt) {
-	const PylithScalar* values = &_data->vertexFields[i][ipt*fiberDim];
-	section->updatePoint(*v_iter, values);
+      PetscSection section = field.petscSection();
+      Vec          vec     = field.localVector();
+      PetscScalar *a;
+      CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+      err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+      for(PetscInt v = vStart; v < vEnd; ++v) {
+        PetscInt dof, off;
+
+        err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < dof; ++d) {
+          a[off+d] = _data->vertexFields[i][(v-vStart)*dof+d];
+        }
       } // for
-      CPPUNIT_ASSERT_EQUAL(_data->numVertices, ipt);
+      err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+      CPPUNIT_ASSERT_EQUAL(_data->numVertices, vEnd-vStart);
     } // for
   } catch (const ALE::Exception& err) {
     throw std::runtime_error(err.msg());
@@ -158,37 +162,40 @@
   try {
     const int nfields = _data->numCellFields;
 
-    const ALE::Obj<topology::SubMesh::SieveMesh>& sieveSubMesh =
-      _submesh->sieveMesh();
-    CPPUNIT_ASSERT(!sieveSubMesh.isNull());
-    const ALE::Obj<topology::SubMesh::SieveMesh::label_sequence>& cells = 
-      sieveSubMesh->heightStratum(1);
-    assert(!cells.isNull());
-    const topology::SubMesh::SieveMesh::label_sequence::iterator cellsBegin = 
-      cells->begin();
-    const topology::SubMesh::SieveMesh::label_sequence::iterator cellsEnd = 
-      cells->end();
+    DM dmMesh = _submesh->dmMesh();
+    PetscInt       cStart, cEnd;
+    const PetscInt height = 1;
+    PetscErrorCode err;
 
+    CPPUNIT_ASSERT(dmMesh);
+    err = DMComplexGetHeightStratum(dmMesh, height, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
     // Set cell fields
     for (int i=0; i < nfields; ++i) {
       const char* name = _data->cellFieldsInfo[i].name;
       const int fiberDim = _data->cellFieldsInfo[i].fiber_dim;
       fields->add(name, name);
       SubMeshField& field = fields->get(name);
-      field.newSection(topology::FieldBase::CELLS_FIELD, fiberDim, 1);
+      field.newSection(topology::FieldBase::CELLS_FIELD, fiberDim, height);
       field.allocate();
       field.vectorFieldType(_data->cellFieldsInfo[i].field_type);
 
-      const ALE::Obj<topology::SubMesh::RealSection>& section = field.section();
-      CPPUNIT_ASSERT(!section.isNull());
-      int icell = 0;
-      for (topology::SubMesh::SieveMesh::label_sequence::iterator c_iter=cellsBegin;
-	   c_iter != cellsEnd;
-	   ++c_iter, ++icell) {
-	const PylithScalar* values = &_data->cellFields[i][icell*fiberDim];
-	section->updatePoint(*c_iter, values);
+      PetscSection section = field.petscSection();
+      Vec          vec     = field.localVector();
+      PetscScalar *a;
+      CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+      err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+      for(PetscInt c = cStart; c < cEnd; ++c) {
+        PetscInt dof, off;
+
+        err = PetscSectionGetDof(section, c, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(section, c, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < dof; ++d) {
+          a[off+d] = _data->cellFields[i][(c-cStart)*dof+d];
+        }
       } // for
-      CPPUNIT_ASSERT_EQUAL(_data->numCells, icell);
+      err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+      CPPUNIT_ASSERT_EQUAL(_data->numCells, cEnd-cStart);
     } // for
   } catch (const ALE::Exception& err) {
     throw std::runtime_error(err.msg());

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterFaultMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterFaultMesh.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterFaultMesh.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -105,15 +105,13 @@
   try {
     const int nfields = _data->numVertexFields;
 
-    const ALE::Obj<topology::SubMesh::SieveMesh>& sieveFaultMesh = 
-      _faultMesh->sieveMesh();
-    CPPUNIT_ASSERT(!sieveFaultMesh.isNull());
-    const ALE::Obj<topology::SubMesh::SieveMesh::label_sequence>& vertices =
-      sieveFaultMesh->depthStratum(0);
-    CPPUNIT_ASSERT(!vertices.isNull());
-    const topology::SubMesh::SieveMesh::label_sequence::iterator verticesEnd =
-      vertices->end();
+    DM dmMesh = _faultMesh->dmMesh();
+    PetscInt       vStart, vEnd;
+    PetscErrorCode err;
 
+    CPPUNIT_ASSERT(dmMesh);
+    err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
     // Set vertex fields
     for (int i=0; i < nfields; ++i) {
       const char* name = _data->vertexFieldsInfo[i].name;
@@ -124,16 +122,22 @@
       field.allocate();
       field.vectorFieldType(_data->vertexFieldsInfo[i].field_type);
 
-      const ALE::Obj<topology::SubMesh::RealSection>& section = field.section();
-      CPPUNIT_ASSERT(!section.isNull());
-      int ipt = 0;
-      for (topology::SubMesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
-	   v_iter != verticesEnd;
-	   ++v_iter, ++ipt) {
-	const PylithScalar* values = &_data->vertexFields[i][ipt*fiberDim];
-	section->updatePoint(*v_iter, values);
+      PetscSection section = field.petscSection();
+      Vec          vec     = field.localVector();
+      PetscScalar *a;
+      CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+      err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+      for(PetscInt v = vStart; v < vEnd; ++v) {
+        PetscInt dof, off;
+
+        err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < dof; ++d) {
+          a[off+d] = _data->vertexFields[i][(v-vStart)*dof+d];
+        }
       } // for
-      CPPUNIT_ASSERT_EQUAL(_data->numVertices, ipt);
+      err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+      CPPUNIT_ASSERT_EQUAL(_data->numVertices, vEnd-vStart);
     } // for
   } catch (const ALE::Exception& err) {
     throw std::runtime_error(err.msg());
@@ -153,17 +157,14 @@
   try {
     const int nfields = _data->numCellFields;
 
-    const ALE::Obj<topology::SubMesh::SieveMesh>& sieveFaultMesh =
-      _faultMesh->sieveMesh();
-    CPPUNIT_ASSERT(!sieveFaultMesh.isNull());
-    const ALE::Obj<topology::SubMesh::SieveMesh::label_sequence>& cells = 
-      sieveFaultMesh->heightStratum(0);
-    assert(!cells.isNull());
-    const topology::SubMesh::SieveMesh::label_sequence::iterator cellsBegin = 
-      cells->begin();
-    const topology::SubMesh::SieveMesh::label_sequence::iterator cellsEnd = 
-      cells->end();
+    DM dmMesh = _faultMesh->dmMesh();
+    PetscInt       cStart, cEnd;
+    const PetscInt height = 0;
+    PetscErrorCode err;
 
+    CPPUNIT_ASSERT(dmMesh);
+    err = DMComplexGetHeightStratum(dmMesh, height, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
     // Set cell fields
     for (int i=0; i < nfields; ++i) {
       const char* name = _data->cellFieldsInfo[i].name;
@@ -174,16 +175,22 @@
       field.allocate();
       field.vectorFieldType(_data->cellFieldsInfo[i].field_type);
 
-      const ALE::Obj<topology::SubMesh::RealSection>& section = field.section();
-      CPPUNIT_ASSERT(!section.isNull());
-      int icell = 0;
-      for (topology::SubMesh::SieveMesh::label_sequence::iterator c_iter=cellsBegin;
-	   c_iter != cellsEnd;
-	   ++c_iter, ++icell) {
-	const PylithScalar* values = &_data->cellFields[i][icell*fiberDim];
-	section->updatePoint(*c_iter, values);
+      PetscSection section = field.petscSection();
+      Vec          vec     = field.localVector();
+      PetscScalar *a;
+      CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+      err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+      for(PetscInt c = cStart; c < cEnd; ++c) {
+        PetscInt dof, off;
+
+        err = PetscSectionGetDof(section, c, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(section, c, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < dof; ++d) {
+          a[off+d] = _data->cellFields[i][(c-cStart)*dof+d];
+        }
       } // for
-      CPPUNIT_ASSERT_EQUAL(_data->numCells, icell);
+      err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+      CPPUNIT_ASSERT_EQUAL(_data->numCells, cEnd-cStart);
     } // for
   } catch (const ALE::Exception& err) {
     throw std::runtime_error(err.msg());

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterMesh.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterMesh.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -98,14 +98,13 @@
   try {
     const int nfields = _data->numVertexFields;
 
-    const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh = _mesh->sieveMesh();
-    CPPUNIT_ASSERT(!sieveMesh.isNull());
-    const ALE::Obj<topology::Mesh::SieveMesh::label_sequence>& vertices =
-      sieveMesh->depthStratum(0);
-    CPPUNIT_ASSERT(!vertices.isNull());
-    const topology::Mesh::SieveMesh::label_sequence::iterator verticesEnd =
-      vertices->end();
+    DM dmMesh = _mesh->dmMesh();
+    PetscInt       vStart, vEnd;
+    PetscErrorCode err;
 
+    CPPUNIT_ASSERT(dmMesh);
+    err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
     // Set vertex fields
     for (int i=0; i < nfields; ++i) {
       const char* name = _data->vertexFieldsInfo[i].name;
@@ -116,16 +115,22 @@
       field.allocate();
       field.vectorFieldType(_data->vertexFieldsInfo[i].field_type);
 
-      const ALE::Obj<topology::Mesh::RealSection>& section = field.section();
-      CPPUNIT_ASSERT(!section.isNull());
-      int ipt = 0;
-      for (topology::Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
-	   v_iter != verticesEnd;
-	   ++v_iter, ++ipt) {
-	const PylithScalar* values = &_data->vertexFields[i][ipt*fiberDim];
-	section->updatePoint(*v_iter, values);
+      PetscSection section = field.petscSection();
+      Vec          vec     = field.localVector();
+      PetscScalar *a;
+      CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+      err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+      for(PetscInt v = vStart; v < vEnd; ++v) {
+        PetscInt dof, off;
+
+        err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < dof; ++d) {
+          a[off+d] = _data->vertexFields[i][(v-vStart)*dof+d];
+        }
       } // for
-      CPPUNIT_ASSERT_EQUAL(_data->numVertices, ipt);
+      err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+      CPPUNIT_ASSERT_EQUAL(_data->numVertices, vEnd-vStart);
     } // for
   } catch (const ALE::Exception& err) {
     throw std::runtime_error(err.msg());
@@ -147,15 +152,21 @@
   try {
     const int nfields = _data->numCellFields;
 
-    const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh = _mesh->sieveMesh();
-    CPPUNIT_ASSERT(!sieveMesh.isNull());
-    const ALE::Obj<topology::Mesh::SieveMesh::label_sequence>& cells = 
-      (0 == _data->cellsLabel) ? 
-      sieveMesh->depthStratum(1) :
-      sieveMesh->getLabelStratum(_data->cellsLabel, _data->labelId);
-    const topology::Mesh::SieveMesh::label_sequence::iterator cellsEnd = 
-      cells->end();
+    DM              dmMesh = _mesh->dmMesh();
+    IS              cellIS = PETSC_NULL;
+    const PetscInt *cells  = PETSC_NULL;
+    PetscInt        cStart, cEnd, numCells;
+    PetscErrorCode  err;
 
+    CPPUNIT_ASSERT(dmMesh);
+    err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+    numCells = cEnd - cStart;
+    if (_data->cellsLabel) {
+      err = DMComplexGetStratumIS(dmMesh, _data->cellsLabel, _data->labelId, &cellIS);CHECK_PETSC_ERROR(err);
+      err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+      err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+    }
+
     // Set cell fields
     for (int i=0; i < nfields; ++i) {
       const char* name = _data->cellFieldsInfo[i].name;
@@ -166,17 +177,27 @@
       field.allocate();
       field.vectorFieldType(_data->cellFieldsInfo[i].field_type);
 
-      const ALE::Obj<topology::Mesh::RealSection>& section = field.section();
-      CPPUNIT_ASSERT(!section.isNull());
-      int icell = 0;
-      for (topology::Mesh::SieveMesh::label_sequence::iterator c_iter=cells->begin();
-	   c_iter != cellsEnd;
-	   ++c_iter, ++icell) {
-	const PylithScalar* values = &_data->cellFields[i][icell*fiberDim];
-	section->updatePoint(*c_iter, values);
+      PetscSection section = field.petscSection();
+      Vec          vec     = field.localVector();
+      PetscScalar *a;
+      CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+      err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+      for(PetscInt c = 0; c < numCells; ++c) {
+        const PetscInt cell = cells ? cells[c] : c+cStart;
+        PetscInt       dof, off;
+
+        err = PetscSectionGetDof(section, cell, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(section, cell, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < dof; ++d) {
+          a[off+d] = _data->cellFields[i][c*dof+d];
+        }
       } // for
-      CPPUNIT_ASSERT_EQUAL(_data->numCells, icell);
+      err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+      CPPUNIT_ASSERT_EQUAL(_data->numCells, numCells);
     } // for
+    if (_data->cellsLabel) {
+      err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+    }
   } catch (const ALE::Exception& err) {
     throw std::runtime_error(err.msg());
   } catch (...) {

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterPoints.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterPoints.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterPoints.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -98,34 +98,39 @@
   try {
     const int nfields = _data->numVertexFields;
 
-    const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh = _mesh->sieveMesh();
-    CPPUNIT_ASSERT(!sieveMesh.isNull());
-    const ALE::Obj<topology::Mesh::SieveMesh::label_sequence>& vertices =
-      sieveMesh->depthStratum(0);
-    CPPUNIT_ASSERT(!vertices.isNull());
-    const topology::Mesh::SieveMesh::label_sequence::iterator verticesEnd =
-      vertices->end();
+    DM dmMesh = _mesh->dmMesh();
+    PetscInt       vStart, vEnd;
+    PetscErrorCode err;
 
+    CPPUNIT_ASSERT(dmMesh);
+    err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
     // Set vertex fields
     for (int i=0; i < nfields; ++i) {
       const char* name = _data->vertexFieldsInfo[i].name;
       const int fiberDim = _data->vertexFieldsInfo[i].fiber_dim;
       fields->add(name, name);
       MeshField& field = fields->get(name);
-      field.newSection(vertices, fiberDim);
+      field.newSection(topology::FieldBase::VERTICES_FIELD, fiberDim);
       field.allocate();
       field.vectorFieldType(_data->vertexFieldsInfo[i].field_type);
 
-      const ALE::Obj<topology::Mesh::RealSection>& section = field.section();
-      CPPUNIT_ASSERT(!section.isNull());
-      int ipt = 0;
-      for (topology::Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
-	   v_iter != verticesEnd;
-	   ++v_iter, ++ipt) {
-	const PylithScalar* values = &_data->vertexFields[i][ipt*fiberDim];
-	section->updatePoint(*v_iter, values);
+      PetscSection section = field.petscSection();
+      Vec          vec     = field.localVector();
+      PetscScalar *a;
+      CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+      err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+      for(PetscInt v = vStart; v < vEnd; ++v) {
+        PetscInt dof, off;
+
+        err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < dof; ++d) {
+          a[off+d] = _data->vertexFields[i][(v-vStart)*dof+d];
+        }
       } // for
-      CPPUNIT_ASSERT_EQUAL(_data->numVertices, ipt);
+      err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+      CPPUNIT_ASSERT_EQUAL(_data->numVertices, vEnd-vStart);
     } // for
   } catch (const ALE::Exception& err) {
     throw std::runtime_error(err.msg());

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterSubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterSubMesh.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterSubMesh.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -112,15 +112,13 @@
   try {
     const int nfields = _data->numVertexFields;
 
-    const ALE::Obj<topology::SubMesh::SieveMesh>& sieveMesh = 
-      _mesh->sieveMesh();
-    CPPUNIT_ASSERT(!sieveMesh.isNull());
-    const ALE::Obj<topology::SubMesh::SieveMesh::label_sequence>& vertices =
-      sieveMesh->depthStratum(0);
-    CPPUNIT_ASSERT(!vertices.isNull());
-    const topology::SubMesh::SieveMesh::label_sequence::iterator verticesEnd =
-      vertices->end();
+    DM dmMesh = _mesh->dmMesh();
+    PetscInt       vStart, vEnd;
+    PetscErrorCode err;
 
+    CPPUNIT_ASSERT(dmMesh);
+    err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
     // Set vertex fields
     for (int i=0; i < nfields; ++i) {
       const char* name = _data->vertexFieldsInfo[i].name;
@@ -131,16 +129,22 @@
       field.allocate();
       field.vectorFieldType(_data->vertexFieldsInfo[i].field_type);
 
-      const ALE::Obj<topology::SubMesh::RealSection>& section = field.section();
-      CPPUNIT_ASSERT(!section.isNull());
-      int ipt = 0;
-      for (topology::SubMesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
-	   v_iter != verticesEnd;
-	   ++v_iter, ++ipt) {
-	const PylithScalar* values = &_data->vertexFields[i][ipt*fiberDim];
-	section->updatePoint(*v_iter, values);
+      PetscSection section = field.petscSection();
+      Vec          vec     = field.localVector();
+      PetscScalar *a;
+      CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+      err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+      for(PetscInt v = vStart; v < vEnd; ++v) {
+        PetscInt dof, off;
+
+        err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < dof; ++d) {
+          a[off+d] = _data->vertexFields[i][(v-vStart)*dof+d];
+        }
       } // for
-      CPPUNIT_ASSERT_EQUAL(_data->numVertices, ipt);
+      err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+      CPPUNIT_ASSERT_EQUAL(_data->numVertices, vEnd-vStart);
     } // for
   } catch (const ALE::Exception& err) {
     throw std::runtime_error(err.msg());
@@ -162,17 +166,14 @@
   try {
     const int nfields = _data->numCellFields;
 
-    const ALE::Obj<topology::SubMesh::SieveMesh>& sieveSubMesh =
-      _submesh->sieveMesh();
-    CPPUNIT_ASSERT(!sieveSubMesh.isNull());
-    const ALE::Obj<topology::SubMesh::SieveMesh::label_sequence>& cells = 
-      sieveSubMesh->heightStratum(1);
-    assert(!cells.isNull());
-    const topology::SubMesh::SieveMesh::label_sequence::iterator cellsBegin = 
-      cells->begin();
-    const topology::SubMesh::SieveMesh::label_sequence::iterator cellsEnd = 
-      cells->end();
+    DM              dmMesh = _submesh->dmMesh();
+    PetscInt        cStart, cEnd, numCells;
+    PetscErrorCode  err;
 
+    CPPUNIT_ASSERT(dmMesh);
+    err = DMComplexGetHeightStratum(dmMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+    numCells = cEnd - cStart;
+
     // Set cell fields
     for (int i=0; i < nfields; ++i) {
       const char* name = _data->cellFieldsInfo[i].name;
@@ -183,16 +184,23 @@
       field.allocate();
       field.vectorFieldType(_data->cellFieldsInfo[i].field_type);
 
-      const ALE::Obj<topology::SubMesh::RealSection>& section = field.section();
-      CPPUNIT_ASSERT(!section.isNull());
-      int icell = 0;
-      for (topology::SubMesh::SieveMesh::label_sequence::iterator c_iter=cellsBegin;
-	   c_iter != cellsEnd;
-	   ++c_iter, ++icell) {
-	const PylithScalar* values = &_data->cellFields[i][icell*fiberDim];
-	section->updatePoint(*c_iter, values);
+      PetscSection section = field.petscSection();
+      Vec          vec     = field.localVector();
+      PetscScalar *a;
+      CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+      err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+      for(PetscInt c = 0; c < numCells; ++c) {
+        const PetscInt cell = c+cStart;
+        PetscInt       dof, off;
+
+        err = PetscSectionGetDof(section, cell, &dof);CHECK_PETSC_ERROR(err);
+        err = PetscSectionGetOffset(section, cell, &off);CHECK_PETSC_ERROR(err);
+        for(PetscInt d = 0; d < dof; ++d) {
+          a[off+d] = _data->cellFields[i][c*dof+d];
+        }
       } // for
-      CPPUNIT_ASSERT_EQUAL(_data->numCells, icell);
+      err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+      CPPUNIT_ASSERT_EQUAL(_data->numCells, numCells);
     } // for
   } catch (const ALE::Exception& err) {
     throw std::runtime_error(err.msg());

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestOutputManager.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestOutputManager.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestOutputManager.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -193,32 +193,35 @@
   iohandler.read(&mesh);
 
   // Set vertex field
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  const ALE::Obj<SieveMesh::label_sequence>& vertices = 
-    sieveMesh->depthStratum(0);
-  CPPUNIT_ASSERT(!vertices.isNull());
-  const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+  DM dmMesh = mesh.dmMesh();
+  PetscInt       vStart, vEnd;
+  PetscErrorCode err;
 
+  CPPUNIT_ASSERT(dmMesh);
+  err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
   MeshField field(mesh);
-  field.newSection(vertices, fiberDim);
+  field.newSection(topology::FieldBase::VERTICES_FIELD, fiberDim);
   field.allocate();
   field.label(label);
   field.vectorFieldType(fieldType);
   field.scale(scale);
-  const ALE::Obj<RealSection>& section = field.section();
-  CPPUNIT_ASSERT(!section.isNull());
+  PetscSection section = field.petscSection();
+  Vec          vec     = field.localVector();
+  PetscScalar *a;
+  CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+  err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+  for(PetscInt v = vStart; v < vEnd; ++v) {
+    PetscInt dof, off;
 
-  CPPUNIT_ASSERT_EQUAL(nvertices, int(vertices->size()));
-  scalar_array values(nvertices*fiberDim);
-  for (int i=0; i < nvertices*fiberDim; ++i)
-    values[i] = fieldValues[i];
-  values /= scale;
-  int ipt = 0;
-  for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
-       v_iter != verticesEnd;
-       ++v_iter, ++ipt) 
-    section->updatePoint(*v_iter, &values[ipt*fiberDim]);
+    err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+    for(PetscInt d = 0; d < dof; ++d) {
+      a[off+d] = fieldValues[(v-vStart)*dof+d]/scale;
+    }
+  } // for
+  err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+  CPPUNIT_ASSERT_EQUAL(nvertices, vEnd-vStart);
 
   spatialdata::geocoords::CSCart cs;
   const int numTimeSteps = 1;
@@ -279,32 +282,37 @@
   iohandler.read(&mesh);
 
   // Set cell field
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  const ALE::Obj<SieveMesh::label_sequence>& cells = 
-    sieveMesh->heightStratum(0);
-  CPPUNIT_ASSERT(!cells.isNull());
-  const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+  DM              dmMesh = mesh.dmMesh();
+  PetscInt        cStart, cEnd, numCells;
+  PetscErrorCode  err;
 
+  CPPUNIT_ASSERT(dmMesh);
+  err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+  numCells = cEnd - cStart;
+
   MeshField field(mesh);
-  field.newSection(cells, fiberDim);
+  field.newSection(topology::FieldBase::CELLS_FIELD, fiberDim);
   field.allocate();
   field.label(label);
   field.vectorFieldType(fieldType);
   field.scale(scale);
-  const ALE::Obj<RealSection>& section = field.section();
-  CPPUNIT_ASSERT(!section.isNull());
+  PetscSection section = field.petscSection();
+  Vec          vec     = field.localVector();
+  PetscScalar *a;
+  CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+  err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+  for(PetscInt c = 0; c < numCells; ++c) {
+    const PetscInt cell = c+cStart;
+    PetscInt       dof, off;
 
-  CPPUNIT_ASSERT_EQUAL(ncells, int(cells->size()));
-  scalar_array values(ncells*fiberDim);
-  for (int i=0; i < ncells*fiberDim; ++i)
-    values[i] = fieldValues[i];
-  values /= scale;
-  int ipt = 0;
-  for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
-       c_iter != cellsEnd;
-       ++c_iter, ++ipt)
-    section->updatePoint(*c_iter, &values[ipt*fiberDim]);
+    err = PetscSectionGetDof(section, cell, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(section, cell, &off);CHECK_PETSC_ERROR(err);
+    for(PetscInt d = 0; d < dof; ++d) {
+      a[off+d] = fieldValues[c*dof+d]/scale;
+    }
+  } // for
+  err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
+  CPPUNIT_ASSERT_EQUAL(ncells, cEnd-cStart);
 
   spatialdata::geocoords::CSCart cs;
   const int numTimeSteps = 1;

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestVertexFilterVecNorm.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestVertexFilterVecNorm.cc	2012-10-02 00:27:11 UTC (rev 20802)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestVertexFilterVecNorm.cc	2012-10-04 12:17:00 UTC (rev 20803)
@@ -82,45 +82,52 @@
   field.vectorFieldType(fieldType);
   field.label(label.c_str());
 
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
-  const ALE::Obj<SieveMesh::label_sequence>& vertices = 
-    sieveMesh->depthStratum(0);
-  CPPUNIT_ASSERT(!vertices.isNull());
-  const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
-  
-  const ALE::Obj<RealSection>& section = field.section();
-  CPPUNIT_ASSERT(!section.isNull());
-  CPPUNIT_ASSERT_EQUAL(nvertices, int(vertices->size()));
-  int ipt = 0;
-  for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
-       v_iter != verticesEnd;
-       ++v_iter, ++ipt) {
-    const PylithScalar* values = &fieldValues[ipt*fiberDim];
-    section->updatePoint(*v_iter, values);
+  DM dmMesh = mesh.dmMesh();
+  PetscInt       vStart, vEnd;
+  PetscErrorCode err;
+
+  CPPUNIT_ASSERT(dmMesh);
+  err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
+  PetscSection section = field.petscSection();
+  Vec          vec     = field.localVector();
+  PetscScalar *a;
+  CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+  CPPUNIT_ASSERT_EQUAL(nvertices, vEnd-vStart);
+  err = VecGetArray(vec, &a);CHECK_PETSC_ERROR(err);
+  for(PetscInt v = vStart; v < vEnd; ++v) {
+    PetscInt dof, off;
+
+    err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+    for(PetscInt d = 0; d < dof; ++d) {
+      a[off+d] = fieldValues[(v-vStart)*dof+d];
+    }
   } // for
+  err = VecRestoreArray(vec, &a);CHECK_PETSC_ERROR(err);
 
   VertexFilterVecNorm<MeshField> filter;
   const MeshField& fieldF = filter.filter(field);
-  const ALE::Obj<RealSection>& sectionF = fieldF.section();
-  CPPUNIT_ASSERT(!sectionF.isNull());
+  PetscSection sectionF = fieldF.petscSection();
+  Vec          vecF     = fieldF.localVector();
+  CPPUNIT_ASSERT(sectionF);CPPUNIT_ASSERT(vecF);
 
   CPPUNIT_ASSERT_EQUAL(fieldTypeE, fieldF.vectorFieldType());
   CPPUNIT_ASSERT_EQUAL(label, std::string(fieldF.label()));
 
-  ipt = 0;
-  for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
-       v_iter != verticesEnd;
-       ++v_iter, ++ipt) {
-    CPPUNIT_ASSERT_EQUAL(fiberDimE, sectionF->getFiberDimension(*v_iter));
-    const PylithScalar* values = sectionF->restrictPoint(*v_iter);
-    CPPUNIT_ASSERT(0 != values);
-    const PylithScalar tolerance = 1.0e-06;
-    for (int i=0; i < fiberDimE; ++i)
-      CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, 
-				   values[i]/fieldValuesE[ipt*fiberDimE+i],
-				   tolerance);
+  const PylithScalar tolerance = 1.0e-06;
+  err = VecGetArray(vecF, &a);CHECK_PETSC_ERROR(err);
+  for(PetscInt v = vStart; v < vEnd; ++v) {
+    PetscInt dof, off;
+
+    err = PetscSectionGetDof(sectionF, v, &dof);CHECK_PETSC_ERROR(err);
+    err = PetscSectionGetOffset(sectionF, v, &off);CHECK_PETSC_ERROR(err);
+    CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
+    for(PetscInt d = 0; d < dof; ++d) {
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, a[off+d]/fieldValuesE[(v-vStart)*dof+d], tolerance);
+    }
   } // for
+  err = VecRestoreArray(vecF, &a);CHECK_PETSC_ERROR(err);
 } // testFilter
 
 



More information about the CIG-COMMITS mailing list