[cig-commits] r22005 - in short/3D/PyLith/trunk/libsrc/pylith: bc faults feassemble materials

[brad at geodynamics.org]

Author: brad
Date: 2013-05-08 17:00:30 -0700 (Wed, 08 May 2013)
New Revision: 22005

Modified:
   short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
   short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
   short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
   short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.hh
   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/GeometryHex3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.hh
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.icc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
Log:
Code cleanup. Added more assert error checking for getClosure(). Remove copies for minCellWidth.

Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -139,8 +139,6 @@
   topology::Field<topology::SubMesh>& dampingConsts = _parameters->get("damping constants");
   topology::VecVisitorMesh dampingConstsVisitor(dampingConsts);
 
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
   topology::CoordsVisitor coordsVisitor(dmSubMesh);
 
   assert(_normalizer);
@@ -158,7 +156,9 @@
 
   for(PetscInt c = cStart; c < cEnd; ++c) {
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
 
     const PetscInt doff = dampingConstsVisitor.sectionOffset(c);
@@ -263,14 +263,10 @@
   // Use _cellVector for cell residual.
   topology::VecVisitorSubMesh residualVisitor(residual, submeshIS);
 
-  PetscScalar *velocityArray = NULL;
-  PetscInt velocitySize = 0;
   topology::VecVisitorSubMesh velocityVisitor(fields->get("velocity(t)"), submeshIS);
 
   submeshIS.deallocate();
   
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
 topology::CoordsVisitor coordsVisitor(dmSubMesh);
 
   // Get 'surface' cells (1 dimension lower than top-level cells)
@@ -288,7 +284,9 @@
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(geometryEvent);
 #endif
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
 
@@ -301,8 +299,9 @@
     _resetCellVector();
 
     // Restrict input fields to cell
-    velocityVisitor.getClosure(&velocityArray, &velocitySize, c);
-    assert(velocitySize == numBasis*spaceDim);
+    PetscScalar *velocityCell = NULL;
+    PetscInt velocitySize = 0;
+    velocityVisitor.getClosure(&velocityCell, &velocitySize, c);assert(velocityCell);assert(numBasis*spaceDim == velocitySize);
 
     const PetscInt doff = dampingConstsVisitor.sectionOffset(c);
     assert(numQuadPts*spaceDim == dampingConstsVisitor.sectionDof(c));
@@ -327,11 +326,11 @@
           for (int iDim=0; iDim < spaceDim; ++iDim)
             _cellVector[iBasis*spaceDim+iDim] -= 
               dampingConstsArray[doff+iQuad*spaceDim+iDim] *
-              valIJ * velocityArray[jBasis*spaceDim+iDim];
+              valIJ * velocityCell[jBasis*spaceDim+iDim];
         } // for
       } // for
     } // for
-    velocityVisitor.restoreClosure(&velocityArray, &velocitySize, c);
+    velocityVisitor.restoreClosure(&velocityCell, &velocitySize, c);
 
     residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), c, ADD_VALUES);
 
@@ -403,14 +402,10 @@
   // Use _cellVector for cell values.
   topology::VecVisitorSubMesh residualVisitor(residual, submeshIS);
 
-  PetscScalar *velocityArray = NULL;
-  PetscInt velocitySize;
   topology::VecVisitorSubMesh velocityVisitor(fields->get("velocity(t)"), submeshIS);
 
   submeshIS.deallocate();
 
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
   topology::CoordsVisitor coordsVisitor(dmSubMesh);
 
   _logger->eventEnd(setupEvent);
@@ -423,7 +418,9 @@
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(geometryEvent);
 #endif
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
 
@@ -436,8 +433,9 @@
     _resetCellVector();
 
     // Restrict input fields to cell
-    velocityVisitor.getClosure(&velocityArray, &velocitySize, c);
-    assert(velocitySize == numBasis*spaceDim);
+    PetscScalar *velocityCell = NULL;
+    PetscInt velocitySize;
+    velocityVisitor.getClosure(&velocityCell, &velocitySize, c);assert(velocityCell);assert(velocitySize == numBasis*spaceDim);
 
     const PetscInt doff = dampingConstsVisitor.sectionOffset(c);
     assert(numQuadPts*spaceDim == dampingConstsVisitor.sectionDof(c));
@@ -464,10 +462,10 @@
         for (int iDim = 0; iDim < spaceDim; ++iDim)
           _cellVector[iBasis*spaceDim+iDim] -= 
             dampingConstsArray[doff+iQuad*spaceDim+iDim] *
-            valIJ * velocityArray[iBasis*spaceDim+iDim];
+            valIJ * velocityCell[iBasis*spaceDim+iDim];
       } // for
     } // for
-    velocityVisitor.restoreClosure(&velocityArray, &velocitySize, c);
+    velocityVisitor.restoreClosure(&velocityCell, &velocitySize, c);
 
     residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), c, ADD_VALUES);
 
@@ -545,8 +543,6 @@
   // Allocate matrix for cell values.
   _initCellMatrix();
 
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
   topology::CoordsVisitor coordsVisitor(dmSubMesh);
 
   _logger->eventEnd(setupEvent);
@@ -559,7 +555,9 @@
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(geometryEvent);
 #endif
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
 
@@ -678,8 +676,6 @@
 
   submeshIS.deallocate();
   
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
   topology::CoordsVisitor coordsVisitor(dmSubMesh);
 
   _logger->eventEnd(setupEvent);
@@ -692,7 +688,9 @@
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(geometryEvent);
 #endif
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -121,13 +121,13 @@
   topology::VecVisitorSubMesh residualVisitor(residual, submeshIS);
   submeshIS.deallocate();
 
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
   topology::CoordsVisitor coordsVisitor(dmSubMesh);
 
   // Loop over faces and integrate contribution from each face
   for(PetscInt c = cStart; c < cEnd; ++c) {
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
 
@@ -417,8 +417,6 @@
   PetscScalar* valueArray = valueVisitor.localArray();
 
   // Get coordinates
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
   topology::CoordsVisitor coordsVisitor(dmSubMesh);
 
   const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
@@ -433,7 +431,9 @@
   // Loop over cells in boundary mesh and perform queries.
   for(PetscInt c = cStart; c < cEnd; ++c) {
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
 
@@ -507,8 +507,6 @@
   scalar_array orientation(orientationSize);
 
   // Get coordinates.
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
   topology::CoordsVisitor coordsVisitor(dmSubMesh);
 
   // Get sections
@@ -530,7 +528,9 @@
   // rotate corresponding traction vector from local to global coordinates.
   for(PetscInt c = cStart; c < cEnd; ++c) {
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
 
     for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts; ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -1818,12 +1818,8 @@
 
   // Get sections
   topology::CoordsVisitor coordsVisitor(faultDMMesh);
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   topology::VecVisitorMesh dLagrangeVisitor(_fields->get("sensitivity dLagrange"));
-  PetscScalar* dLagrangeCell = NULL;
-  PetscInt dLagrangeSize = 0;
 
   scalar_array residualCell(numBasis*spaceDim);
   topology::Field<topology::SubMesh>& residual = _fields->get("sensitivity residual");
@@ -1833,12 +1829,16 @@
   // Loop over cells
   for(PetscInt c = cStart; c < cEnd; ++c) {
     // Compute geometry
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
 
     // Restrict input fields to cell
-    dLagrangeVisitor.getClosure(&dLagrangeCell, &dLagrangeSize, c);
+    PetscScalar* dLagrangeCell = NULL;
+    PetscInt dLagrangeSize = 0;
+    dLagrangeVisitor.getClosure(&dLagrangeCell, &dLagrangeSize, c);assert(dLagrangeCell);assert(numBasis*spaceDim == dLagrangeSize);
 
     // Get cell geometry information that depends on cell
     const scalar_array& basis = _quadrature->basis();

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -1310,8 +1310,6 @@
 
   // Get section containing coordinates of vertices
   topology::CoordsVisitor coordsVisitor(faultDMMesh);
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   // Loop over cohesive cells, computing orientation weighted by
   // jacobian at constraint vertices
@@ -1322,7 +1320,9 @@
 
     // Get orientations at fault cell's vertices.
 
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
 
     PetscErrorCode err = DMPlexGetTransitiveClosure(faultDMMesh, c, PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
 
@@ -1562,15 +1562,15 @@
   scalar_array areaCell(numBasis);
 
   topology::CoordsVisitor coordsVisitor(faultDMMesh);
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   // Loop over cells in fault mesh, compute area
   for(PetscInt c = cStart; c < cEnd; ++c) {
     areaCell = 0.0;
 
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, c);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, c);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -160,10 +160,15 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
   virtual
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const = 0;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const = 0;
 
   /** Compute orientation of cell at location.
    *

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -207,23 +207,15 @@
 
   // Setup field visitors.
   topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
-  PetscScalar* accCell = NULL;
-  PetscInt accSize = 0;
 
   topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
-  PetscScalar* velCell = NULL;
-  PetscInt velSize = 0;
 
   scalar_array dispAdjCell(numBasis*spaceDim);
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-  PetscScalar* dispCell = NULL;
-  PetscInt dispSize = 0;
   
   topology::VecVisitorMesh residualVisitor(residual);
 
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   assert(_normalizer);
   const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -247,7 +239,9 @@
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(geometryEvent);
 #endif
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 
@@ -268,12 +262,18 @@
     _resetCellVector();
 
     // Restrict input fields to cell
-    accVisitor.getClosure(&accCell, &accSize, cell);
-    velVisitor.getClosure(&velCell, &velSize, cell);
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);
-    assert(velSize == accSize);
-    assert(dispSize == accSize);
+    PetscScalar* accCell = NULL;
+    PetscInt accSize = 0;
+    accVisitor.getClosure(&accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
 
+    PetscScalar* velCell = NULL;
+    PetscInt velSize = 0;
+    velVisitor.getClosure(&velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
+
+    PetscScalar* dispCell = NULL;
+    PetscInt dispSize = 0;
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(restrictEvent);
     _logger->eventBegin(computeEvent);
@@ -450,8 +450,6 @@
   PetscInt jacobianSize = 0;
 
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   _logger->eventEnd(setupEvent);
 #if !defined(DETAILED_EVENT_LOGGING)
@@ -465,7 +463,9 @@
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(geometryEvent);
 #endif
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -200,24 +200,16 @@
 
   // Setup field visitors.
   topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
-  PetscScalar* accCell = NULL;
-  PetscInt accSize = 0;
 
   topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
-  PetscScalar* velCell = NULL;
-  PetscInt velSize = 0;
 
   scalar_array dispAdjCell(numBasis*spaceDim);
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-  PetscScalar* dispCell = NULL;
-  PetscInt dispSize = 0;
   
   topology::VecVisitorMesh residualVisitor(residual);
 
   scalar_array coordinatesCell(numBasis*spaceDim);
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   assert(_normalizer);
   const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -235,8 +227,11 @@
   // Loop over cells
   for(PetscInt c = 0; c < numCells; ++c) {
     const PetscInt cell = cells[c];
+
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
 
     // Get state variables for cell.
@@ -246,12 +241,18 @@
     _resetCellVector();
 
     // Restrict input fields to cell
-    accVisitor.getClosure(&accCell, &accSize, cell);
-    velVisitor.getClosure(&velCell, &velSize, cell);
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);
-    assert(velSize == accSize);
-    assert(dispSize == accSize);
+    PetscScalar* accCell = NULL;
+    PetscInt accSize = 0;
+    accVisitor.getClosure(&accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
 
+    PetscScalar* velCell = NULL;
+    PetscInt velSize = 0;
+    velVisitor.getClosure(&velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
+
+    PetscScalar* dispCell = NULL;
+    PetscInt dispSize = 0;
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+
     // Get cell geometry information that depends on cell
     const scalar_array& basis = _quadrature->basis();
     const scalar_array& basisDeriv = _quadrature->basisDeriv();
@@ -402,8 +403,6 @@
   topology::VecVisitorMesh jacobianVisitor(*jacobian);
 
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   _logger->eventEnd(setupEvent);
   _logger->eventBegin(computeEvent);
@@ -412,7 +411,9 @@
   for(PetscInt c = 0; c < numCells; ++c) {
     const PetscInt cell = cells[c];
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -192,23 +192,15 @@
 
   // Setup field visitors.
   topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
-  PetscScalar* accCell = NULL;
-  PetscInt accSize = 0;
 
   topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
-  PetscScalar* velCell = NULL;
-  PetscInt velSize = 0;
 
   scalar_array dispAdjCell(numBasis*spaceDim);
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-  PetscScalar* dispCell = NULL;
-  PetscInt dispSize = 0;
   
   topology::VecVisitorMesh residualVisitor(residual);
 
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   assert(_normalizer);
   const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -233,11 +225,17 @@
 #endif
 
     // Restrict input fields to cell
-    accVisitor.getClosure(&accCell, &accSize, cell);
-    velVisitor.getClosure(&velCell, &velSize, cell);
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);
-    assert(velSize == accSize);
-    assert(dispSize == accSize);
+    PetscScalar* accCell = NULL;
+    PetscInt accSize = 0;
+    accVisitor.getClosure(&accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
+    
+    PetscScalar* velCell = NULL;
+    PetscInt velSize = 0;
+    velVisitor.getClosure(&velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
+    
+    PetscScalar* dispCell = NULL;
+    PetscInt dispSize = 0;
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
 
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(restrictEvent);
@@ -245,7 +243,9 @@
 #endif
 
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     const PylithScalar volume = _volume(coordsCell, coordsSize);assert(volume > 0.0);
 
 #if defined(DETAILED_EVENT_LOGGING)
@@ -491,8 +491,6 @@
   PetscInt jacobianSize = 0;
 
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   _logger->eventEnd(setupEvent);
 #if !defined(DETAILED_EVENT_LOGGING)
@@ -505,7 +503,9 @@
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(geometryEvent);
 #endif
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     const PylithScalar volume = _volume(coordsCell, coordsSize);assert(volume > 0.0);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -191,23 +191,15 @@
 
   // Setup field visitors.
   topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
-  PetscScalar* accCell = NULL;
-  PetscInt accSize = 0;
 
   topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
-  PetscScalar* velCell = NULL;
-  PetscInt velSize = 0;
 
   scalar_array dispAdjCell(numBasis*spaceDim);
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-  PetscScalar* dispCell = NULL;
-  PetscInt dispSize = 0;
   
   topology::VecVisitorMesh residualVisitor(residual);
 
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   assert(_normalizer);
   const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -232,19 +224,27 @@
 #endif
 
     // Restrict input fields to cell
-    accVisitor.getClosure(&accCell, &accSize, cell);
-    velVisitor.getClosure(&velCell, &velSize, cell);
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);
-    assert(velSize == accSize);
-    assert(dispSize == accSize);
+    PetscScalar* accCell = NULL;
+    PetscInt accSize = 0;
+    accVisitor.getClosure(&accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
 
+    PetscScalar* velCell = NULL;
+    PetscInt velSize = 0;
+    velVisitor.getClosure(&velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
+
+    PetscScalar* dispCell = NULL;
+    PetscInt dispSize = 0;
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(restrictEvent);
     _logger->eventBegin(geometryEvent);
 #endif
 
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     const PylithScalar area = _area(coordsCell, coordsSize);assert(area > 0.0);
 
 #if defined(DETAILED_EVENT_LOGGING)
@@ -453,8 +453,6 @@
   PetscInt jacobianSize = 0;
 
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   _logger->eventEnd(setupEvent);
 #if !defined(DETAILED_EVENT_LOGGING)
@@ -467,7 +465,9 @@
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(geometryEvent);
 #endif
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     const PylithScalar area = _area(coordsCell, coordsSize);assert(area > 0.0);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -174,18 +174,9 @@
 
   // Setup field visitors.
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-  PetscScalar* dispCell = NULL;
-  PetscInt dispSize = 0;
-
   topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
-  PetscScalar* dispIncrCell = NULL;
-  PetscInt dispIncrSize = 0;
-
   topology::VecVisitorMesh residualVisitor(residual);
-
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   assert(_normalizer);
   const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -201,7 +192,9 @@
 #if defined(PRECOMPUTE_GEOMETRY)
     _quadrature->retrieveGeometry(cell);
 #else
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 #endif
@@ -213,10 +206,14 @@
     _resetCellVector();
 
     // Restrict input fields to cell
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);
-    dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);
-    assert(dispSize == dispIncrSize);
+    PetscScalar* dispCell = NULL;
+    PetscInt dispSize = 0;
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
 
+    PetscScalar* dispIncrCell = NULL;
+    PetscInt dispIncrSize = 0;
+    dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);assert(numBasis*spaceDim == dispIncrSize);
+
     // Get cell geometry information that depends on cell
     const scalar_array& basis = _quadrature->basis();
     const scalar_array& basisDeriv = _quadrature->basisDeriv();
@@ -358,16 +355,8 @@
 
   // Setup field visitors.
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-  PetscScalar* dispCell = NULL;
-  PetscInt dispSize = 0;
-
   topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
-  PetscScalar* dispIncrCell = NULL;
-  PetscInt dispIncrSize = 0;
-
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   // Get sparse matrix
   const PetscMat jacobianMat = jacobian->matrix();assert(jacobianMat);
@@ -387,7 +376,9 @@
 #if defined(PRECOMPUTE_GEOMETRY)
     _quadrature->retrieveGeometry(cell);
 #else
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 #endif
@@ -399,10 +390,14 @@
     _resetCellMatrix();
 
     // Restrict input fields to cell
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);
-    dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);
-    assert(dispSize == dispIncrSize);
+    PetscScalar* dispCell = NULL;
+    PetscInt dispSize = 0;
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
 
+    PetscScalar* dispIncrCell = NULL;
+    PetscInt dispIncrSize = 0;
+    dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);assert(numBasis*spaceDim == dispIncrSize);
+
     // Get cell geometry information that depends on cell
     const scalar_array& basis = _quadrature->basis();
     const scalar_array& basisDeriv = _quadrature->basisDeriv();

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -178,11 +178,8 @@
 
   // Setup field visitors.
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-
   topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
-
   topology::VecVisitorMesh residualVisitor(residual);
-
   topology::CoordsVisitor coordsVisitor(dmMesh);
 
   assert(_normalizer);
@@ -213,12 +210,11 @@
     // Restrict input fields to cell
     PetscScalar* dispCell = NULL;
     PetscInt dispSize = 0;
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
 
     PetscScalar* dispIncrCell = NULL;
     PetscInt dispIncrSize = 0;
-    dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);
-    assert(dispSize == dispIncrSize);
+    dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);assert(numBasis*spaceDim == dispIncrSize);
 
     // Get cell geometry information that depends on cell
     const scalar_array& basis = _quadrature->basis();
@@ -366,17 +362,8 @@
 
   // Setup field visitors.
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-  PetscScalar* dispCell = NULL;
-  PetscInt dispSize = 0;
-
   topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
-  PetscScalar* dispIncrCell = NULL;
-  PetscInt dispIncrSize = 0;
-
-  scalar_array coordinatesCell(numBasis*spaceDim);
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   // Get sparse matrix
   const PetscMat jacobianMat = jacobian->matrix();assert(jacobianMat);
@@ -393,10 +380,11 @@
   for(PetscInt c = 0; c < numCells; ++c) {
     const PetscInt cell = cells[c];
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
 
-
     // Get state variables for cell.
     _material->retrievePropsAndVars(cell);
 
@@ -404,10 +392,14 @@
     _resetCellMatrix();
 
     // Restrict input fields to cell
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);
-    dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);
-    assert(dispSize == dispIncrSize);
+    PetscScalar* dispCell = NULL;
+    PetscInt dispSize = 0;
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
 
+    PetscScalar* dispIncrCell = NULL;
+    PetscInt dispIncrSize = 0;
+    dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);assert(numBasis*spaceDim == dispIncrSize);
+
     // Get cell geometry information that depends on cell
     const scalar_array& basis = _quadrature->basis();
     const scalar_array& basisDeriv = _quadrature->basisDeriv();

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -378,12 +378,17 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryHex3D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryHex3D::minCellWidth(const PylithScalar* coordinatesCell,
+						const int numVertices,
+						const int spaceDim) const
 { // minCellWidth
   const int numCorners = 8;
-  const int spaceDim = 3;
-  assert(numCorners*spaceDim == coordinatesCell.size());
+  const int dim = 3;
 
+  assert(coordinatesCell);
+  assert(numCorners == numVertices);
+  assert(dim == spaceDim);
+
   const int numEdges = 12;
   const int edges[numEdges][2] = {
     {0, 1}, {1, 2}, {2, 3}, {3, 0},
@@ -395,12 +400,12 @@
   for (int iedge=0; iedge < numEdges; ++iedge) {
     const int iA = edges[iedge][0];
     const int iB = edges[iedge][1];
-    const PylithScalar xA = coordinatesCell[spaceDim*iA  ];
-    const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
-    const PylithScalar zA = coordinatesCell[spaceDim*iA+2];
-    const PylithScalar xB = coordinatesCell[spaceDim*iB  ];
-    const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
-    const PylithScalar zB = coordinatesCell[spaceDim*iB+2];
+    const PylithScalar xA = coordinatesCell[dim*iA  ];
+    const PylithScalar yA = coordinatesCell[dim*iA+1];
+    const PylithScalar zA = coordinatesCell[dim*iA+2];
+    const PylithScalar xB = coordinatesCell[dim*iB  ];
+    const PylithScalar yB = coordinatesCell[dim*iB+1];
+    const PylithScalar zB = coordinatesCell[dim*iB+2];
     
     const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2) + pow(zB-zA,2));
     if (edgeLen < minWidth) {

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -123,9 +123,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -151,13 +151,18 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryLine1D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryLine1D::minCellWidth(const PylithScalar* coordinatesCell,
+						 const int numVertices,
+						 const int spaceDim) const
 { // minCellWidth
   const int numCorners = 2;
-  const int spaceDim = 1;
-  assert(2*spaceDim == coordinatesCell.size() ||
-	 3*spaceDim == coordinatesCell.size()); // :KLUDGE: allow quadratic
+  const int dim = 1;
 
+  assert(coordinatesCell);
+  assert(numCorners == numVertices || // linear
+	 numCorners+1 == numVertices); // quadratic
+  assert(dim == spaceDim);
+
   const PylithScalar xA = coordinatesCell[0];
   const PylithScalar xB = coordinatesCell[1];
     

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -107,9 +107,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -172,12 +172,16 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryLine2D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryLine2D::minCellWidth(const PylithScalar* coordinatesCell,
+						 const int numVertices,
+						 const int spaceDim) const
 { // minCellWidth
   const int numCorners = 2;
-  const int spaceDim = 2;
-  assert(numCorners*spaceDim == coordinatesCell.size());
+  const int dim = 2;
 
+  assert(numCorners == numVertices);
+  assert(dim == spaceDim);
+
   const PylithScalar xA = coordinatesCell[0];
   const PylithScalar yA = coordinatesCell[1];
   const PylithScalar xB = coordinatesCell[2];

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -107,9 +107,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -183,12 +183,16 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryLine3D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryLine3D::minCellWidth(const PylithScalar* coordinatesCell,
+						 const int numVertices,
+						 const int spaceDim) const
 { // minCellWidth
   const int numCorners = 2;
-  const int spaceDim = 3;
-  assert(numCorners*spaceDim == coordinatesCell.size());
+  const int dim = 3;
 
+  assert(numCorners == numVertices);
+  assert(dim == spaceDim);
+
   const PylithScalar xA = coordinatesCell[0];
   const PylithScalar yA = coordinatesCell[1];
   const PylithScalar zA = coordinatesCell[2];

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -107,9 +107,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -122,7 +122,9 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryPoint1D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryPoint1D::minCellWidth(const PylithScalar* coordinatesCell,
+						  const int numVertices,
+						  const int spaceDim) const
 { // minCellWidth
   return PYLITH_MAXSCALAR;
 } // minCellWidth

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -107,9 +107,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -123,7 +123,9 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryPoint2D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryPoint2D::minCellWidth(const PylithScalar* coordinatesCell,
+						  const int numVertices,
+						  const int spaceDim) const
 { // minCellWidth
   return PYLITH_MAXSCALAR;
 } // minCellWidth

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -107,9 +107,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -123,7 +123,9 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryPoint3D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryPoint3D::minCellWidth(const PylithScalar* coordinatesCell,
+						  const int numVertices,
+						  const int spaceDim) const
 { // minCellWidth
   return PYLITH_MAXSCALAR;
 } // minCellWidth

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -107,9 +107,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -227,12 +227,16 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryQuad2D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryQuad2D::minCellWidth(const PylithScalar* coordinatesCell,
+						 const int numVertices,
+						 const int spaceDim) const
 { // minCellWidth
   const int numCorners = 4;
-  const int spaceDim = 2;
-  assert(numCorners*spaceDim == coordinatesCell.size());
+  const int dim = 2;
 
+  assert(numCorners == numVertices);
+  assert(dim == spaceDim);
+
   const int numEdges = 4;
   const int edges[numEdges][2] = {
     {0, 1}, {1, 2}, {2, 3}, {3, 0},
@@ -242,10 +246,10 @@
   for (int iedge=0; iedge < numEdges; ++iedge) {
     const int iA = edges[iedge][0];
     const int iB = edges[iedge][1];
-    const PylithScalar xA = coordinatesCell[spaceDim*iA  ];
-    const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
-    const PylithScalar xB = coordinatesCell[spaceDim*iB  ];
-    const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
+    const PylithScalar xA = coordinatesCell[dim*iA  ];
+    const PylithScalar yA = coordinatesCell[dim*iA+1];
+    const PylithScalar xB = coordinatesCell[dim*iB  ];
+    const PylithScalar yB = coordinatesCell[dim*iB+1];
     
     const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2));
     if (edgeLen < minWidth) {

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -125,9 +125,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -272,12 +272,16 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryQuad3D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryQuad3D::minCellWidth(const PylithScalar* coordinatesCell,
+						 const int numVertices,
+						 const int spaceDim) const
 { // minCellWidth
   const int numCorners = 4;
-  const int spaceDim = 3;
-  assert(numCorners*spaceDim == coordinatesCell.size());
+  const int dim = 3;
 
+  assert(numCorners == numVertices);
+  assert(dim == spaceDim);
+
   const int numEdges = 4;
   const int edges[numEdges][2] = {
     {0, 1}, {1, 2}, {2, 3}, {3, 0},
@@ -287,12 +291,12 @@
   for (int iedge=0; iedge < numEdges; ++iedge) {
     const int iA = edges[iedge][0];
     const int iB = edges[iedge][1];
-    const PylithScalar xA = coordinatesCell[spaceDim*iA  ];
-    const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
-    const PylithScalar zA = coordinatesCell[spaceDim*iA+2];
-    const PylithScalar xB = coordinatesCell[spaceDim*iB  ];
-    const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
-    const PylithScalar zB = coordinatesCell[spaceDim*iB+2];
+    const PylithScalar xA = coordinatesCell[dim*iA  ];
+    const PylithScalar yA = coordinatesCell[dim*iA+1];
+    const PylithScalar zA = coordinatesCell[dim*iA+2];
+    const PylithScalar xB = coordinatesCell[dim*iB  ];
+    const PylithScalar yB = coordinatesCell[dim*iB+1];
+    const PylithScalar zB = coordinatesCell[dim*iB+2];
     
     const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2) + pow(zB-zA,2));
     if (edgeLen < minWidth) {

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -123,9 +123,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -253,13 +253,17 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryTet3D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryTet3D::minCellWidth(const PylithScalar* coordinatesCell,
+						const int numVertices,
+						const int spaceDim) const
 { // minCellWidth
   const int numCorners = 4;
-  const int spaceDim = 3;
-  assert(4*spaceDim == coordinatesCell.size() ||
-	 10*spaceDim == coordinatesCell.size()); // :KLUDGE: allow quadratic
+  const int dim = 3;
 
+  assert(numCorners == numVertices || // linear
+	 numCorners+6 == numVertices); // quadratic
+  assert(dim == spaceDim);
+
   const int numEdges = 6;
   const int edges[numEdges][2] = {
     {0, 1}, {1, 2}, {2, 0},
@@ -270,12 +274,12 @@
   for (int iedge=0; iedge < numEdges; ++iedge) {
     const int iA = edges[iedge][0];
     const int iB = edges[iedge][1];
-    const PylithScalar xA = coordinatesCell[spaceDim*iA  ];
-    const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
-    const PylithScalar zA = coordinatesCell[spaceDim*iA+2];
-    const PylithScalar xB = coordinatesCell[spaceDim*iB  ];
-    const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
-    const PylithScalar zB = coordinatesCell[spaceDim*iB+2];
+    const PylithScalar xA = coordinatesCell[dim*iA  ];
+    const PylithScalar yA = coordinatesCell[dim*iA+1];
+    const PylithScalar zA = coordinatesCell[dim*iA+2];
+    const PylithScalar xB = coordinatesCell[dim*iB  ];
+    const PylithScalar yB = coordinatesCell[dim*iB+1];
+    const PylithScalar zB = coordinatesCell[dim*iB+2];
     
     const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2) + pow(zB-zA,2));
     if (edgeLen < minWidth) {
@@ -286,11 +290,11 @@
   PetscLogFlops(numEdges*9);
 
   // Radius of inscribed sphere
-  const PylithScalar v = volume(coordinatesCell);
-  const PylithScalar a = faceArea(coordinatesCell, 0) +
-    faceArea(coordinatesCell, 1) +
-    faceArea(coordinatesCell, 2) +
-    faceArea(coordinatesCell, 3);
+  const PylithScalar v = volume(coordinatesCell, numVertices, spaceDim);
+  const PylithScalar a = faceArea(coordinatesCell, numVertices, spaceDim, 0) +
+    faceArea(coordinatesCell, numVertices, spaceDim, 1) +
+    faceArea(coordinatesCell, numVertices, spaceDim, 2) +
+    faceArea(coordinatesCell, numVertices, spaceDim, 3);
     
   const PylithScalar r = 3.0 * v / a;
   const PylithScalar rwidth = 6.38*r; // based on empirical tests
@@ -306,10 +310,16 @@
 // ----------------------------------------------------------------------
 // Compute cell volume.
 PylithScalar
-pylith::feassemble::GeometryTet3D::volume(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryTet3D::volume(const PylithScalar* coordinatesCell,
+					  const int numVertices,
+					  const int spaceDim) const
 { // volume
-  assert(12 == coordinatesCell.size() ||
-	 30 == coordinatesCell.size()); // :KLUDGE: allow quadratic
+  const int numCorners = 4;
+  const int dim = 3;
+
+  assert(numCorners == numVertices || // linear
+	 numCorners+6 == numVertices); // quadratic
+  assert(dim == spaceDim);
   
   const PylithScalar x0 = coordinatesCell[0];
   const PylithScalar y0 = coordinatesCell[1];
@@ -342,12 +352,18 @@
 // ----------------------------------------------------------------------
 // Compute area of face.
 PylithScalar
-pylith::feassemble::GeometryTet3D::faceArea(const scalar_array& coordinatesCell,
-	 const int face) const
+pylith::feassemble::GeometryTet3D::faceArea(const PylithScalar* coordinatesCell,
+					    const int numVertices,
+					    const int spaceDim,
+					    const int face) const
 { // faceArea
-  assert(12 == coordinatesCell.size() ||
-	 30 == coordinatesCell.size()); // :KLUDGE: allow quadratic
+  const int numCorners = 4;
+  const int dim = 3;
 
+  assert(numCorners == numVertices || // linear
+	 numCorners+6 == numVertices); // quadratic
+  assert(dim == spaceDim);
+
   const PylithScalar x0 = coordinatesCell[0];
   const PylithScalar y0 = coordinatesCell[1];
   const PylithScalar z0 = coordinatesCell[2];

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -118,24 +118,37 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   /** Compute cell volume.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
    * @returns Volume of cell.
    */
-  PylithScalar volume(const scalar_array& coordinatesCell) const;
+  PylithScalar volume(const PylithScalar* coordinatesCell,
+		      const int numVertices,
+		      const int spaceDim) const;
 
   /** Compute area of face.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
    * @param face Index of vertex across from face.
    * @returns Area of cell face.
    */
-  PylithScalar faceArea(const scalar_array& coordinatesCell,
+  PylithScalar faceArea(const PylithScalar* coordinatesCell,
+			const int numVertices,
+			const int spaceDim,
 			const int face) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -195,13 +195,17 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryTri2D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryTri2D::minCellWidth(const PylithScalar* coordinatesCell,
+						const int numVertices,
+						const int spaceDim) const
 { // minCellWidth
   const int numCorners = 3;
-  const int spaceDim = 2;
-  assert(3*spaceDim == coordinatesCell.size() ||
-	 6*spaceDim == coordinatesCell.size()); // :KLUDGE: allow quadratic
+  const int dim = 2;
 
+  assert(numCorners == numVertices || // linear
+	 numCorners+3 == numVertices); // quadratic
+  assert(dim == spaceDim);
+
   const int numEdges = 3;
   const int edges[numEdges][2] = {
     {0, 1}, {1, 2}, {2, 0},
@@ -211,10 +215,10 @@
   for (int iedge=0; iedge < numEdges; ++iedge) {
     const int iA = edges[iedge][0];
     const int iB = edges[iedge][1];
-    const PylithScalar xA = coordinatesCell[spaceDim*iA  ];
-    const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
-    const PylithScalar xB = coordinatesCell[spaceDim*iB  ];
-    const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
+    const PylithScalar xA = coordinatesCell[dim*iA  ];
+    const PylithScalar yA = coordinatesCell[dim*iA+1];
+    const PylithScalar xB = coordinatesCell[dim*iB  ];
+    const PylithScalar yB = coordinatesCell[dim*iB+1];
     
     const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2));
     if (edgeLen < minWidth) {

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -124,9 +124,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -231,12 +231,16 @@
 // ----------------------------------------------------------------------
 // Compute minimum width across cell.
 PylithScalar
-pylith::feassemble::GeometryTri3D::minCellWidth(const scalar_array& coordinatesCell) const
+pylith::feassemble::GeometryTri3D::minCellWidth(const PylithScalar* coordinatesCell,
+						const int numVertices,
+						const int spaceDim) const
 { // minCellWidth
   const int numCorners = 2;
-  const int spaceDim = 3;
-  assert(numCorners*spaceDim == coordinatesCell.size());
+  const int dim = 3;
 
+  assert(numCorners == numVertices);
+  assert(dim == spaceDim);
+
   const int numEdges = 3;
   const int edges[numEdges][2] = {
     {0, 1}, {1, 2}, {2, 0},
@@ -246,12 +250,12 @@
   for (int iedge=0; iedge < numEdges; ++iedge) {
     const int iA = edges[iedge][0];
     const int iB = edges[iedge][1];
-    const PylithScalar xA = coordinatesCell[spaceDim*iA  ];
-    const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
-    const PylithScalar zA = coordinatesCell[spaceDim*iA+2];
-    const PylithScalar xB = coordinatesCell[spaceDim*iB  ];
-    const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
-    const PylithScalar zB = coordinatesCell[spaceDim*iB+2];
+    const PylithScalar xA = coordinatesCell[dim*iA  ];
+    const PylithScalar yA = coordinatesCell[dim*iA+1];
+    const PylithScalar zA = coordinatesCell[dim*iA+2];
+    const PylithScalar xB = coordinatesCell[dim*iB  ];
+    const PylithScalar yB = coordinatesCell[dim*iB+1];
+    const PylithScalar zB = coordinatesCell[dim*iB+2];
     
     const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2) + pow(zB-zA,2));
     if (edgeLen < minWidth) {

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -124,9 +124,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -515,27 +515,26 @@
 
   // Setup field visitors.
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-  PetscScalar* dispCell = NULL;
-  PetscInt dispSize = 0;
 
   topology::VecVisitorMesh fieldVisitor(*field);
   PetscScalar* fieldArray = fieldVisitor.localArray();
 
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar *coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   // Loop over cells
   for(PetscInt c = 0; c < numCells; ++c) {
     const PetscInt cell = cells[c];
     // Retrieve geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
     coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 
     // Get cell geometry information that depends on cell
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);
-    assert(numBasis*spaceDim == dispSize);
+    PetscScalar* dispCell = NULL;
+    PetscInt dispSize = 0;
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
     const scalar_array& basisDeriv = _quadrature->basisDeriv();
     
     // Compute strains

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -111,24 +111,23 @@
 
   scalar_array dispTCell(numBasis*spaceDim);
   topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
-  PetscScalar* dispCell = NULL;
-  PetscInt dispSize = 0;
 
   scalar_array coordinatesCell(numBasis*spaceDim);
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar* coordsCell = NULL;
-  PetscInt coordsSize = 0;
 
   // Loop over cells
   for (PetscInt c = 0; c < numCells; ++c) {
     const PetscInt cell = cells[c];
     // Retrieve geometry information for current cell
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
     const scalar_array& basisDeriv = _quadrature->basisDeriv();
 
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);
-    assert(numBasis*spaceDim == dispSize);
+    PetscScalar* dispCell = NULL;
+    PetscInt dispSize = 0;
+    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
   
     // Compute deformation tensor.
     _calcDeformation(&deformCell, basisDeriv, coordsCell, dispCell, numBasis, numQuadPts, spaceDim);

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.hh	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.hh	2013-05-09 00:00:30 UTC (rev 22005)
@@ -190,9 +190,14 @@
   /** Compute minimum width across cell.
    *
    * @param coordinatesCell Coordinates of vertices in cell.
+   * @param numVertices Number of vertices in cell.
+   * @param spaceDim Coordinate dimension.
+   *
    * @returns Minimum width across cell.
    */
-  PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+  PylithScalar minCellWidth(const PylithScalar* coordinatesCell,
+			    const int numVertices,
+			    const int spaceDim) const;
 
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.icc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.icc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -93,9 +93,11 @@
 // Compute minimum width across cell.
 inline
 PylithScalar
-pylith::feassemble::QuadratureRefCell::minCellWidth(const scalar_array& coordinatesCell) const {
+pylith::feassemble::QuadratureRefCell::minCellWidth(const PylithScalar* coordinatesCell,
+						    const int numVertices,
+						    const int spaceDim) const {
   assert(_geometry);
-  return _geometry->minCellWidth(coordinatesCell);
+  return _geometry->minCellWidth(coordinatesCell, numVertices, spaceDim);
 }
 
 #endif

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -416,8 +416,6 @@
 
   scalar_array coordinatesCell(numBasis*spaceDim);
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar *coordsArray = NULL;
-  PetscInt coordsSize = 0;
 
   PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
   scalar_array dtStableCell(numQuadPts);
@@ -426,12 +424,11 @@
 
     retrievePropsAndVars(cell);
 
-    coordsVisitor.getClosure(&coordsArray, &coordsSize, cell);
-    for(PetscInt i = 0; i < coordsSize; ++i) {
-      coordinatesCell[i] = coordsArray[i];
-    } // for
-    coordsVisitor.restoreClosure(&coordsArray, &coordsSize, cell);
-    const double minCellWidth = quadrature->minCellWidth(coordinatesCell);
+    PetscScalar *coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
+    const PylithScalar minCellWidth = quadrature->minCellWidth(coordsCell, numBasis, spaceDim);
+    coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 
     for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
       const PylithScalar dt = 
@@ -590,7 +587,7 @@
 
 #if !defined(PRECOMPUTE_GEOMETRY)
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar *coordsArray = NULL;
+  PetscScalar *coordsCell = NULL;
   PetscInt coordsSize = 0;
 #endif
 
@@ -654,9 +651,9 @@
 #if defined(PRECOMPUTE_GEOMETRY)
     quadrature->retrieveGeometry(*c_iter);
 #else
-    coordsVisitor.getClosure(&coordsArray, &coordsSize, cell);
-    quadrature->computeGeometry(coordsArray, coordsSize, cell);
-    coordsVisitor.restoreClosure(&coordsArray, &coordsSize, cell);
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
+    quadrature->computeGeometry(coordsCell, coordsSize, cell);
+    coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 #endif
 
     // Dimensionalize coordinates for querying
@@ -728,7 +725,7 @@
 
 #if !defined(PRECOMPUTE_GEOMETRY)
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar* coordsArray = NULL;
+  PetscScalar* coordsCell = NULL;
   PetscInt coordsSize = 0;
 #endif
 
@@ -792,9 +789,9 @@
 #if defined(PRECOMPUTE_GEOMETRY)
     quadrature->retrieveGeometry(*c_iter);
 #else
-    coordsVisitor.getClosure(&coordsArray, &coordsSize, cell);
-    quadrature->computeGeometry(coordsArray, coordsSize, cell);
-    coordsVisitor.restoreClosure(&coordsArray, &coordsSize, cell);
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
+    quadrature->computeGeometry(coordsCell, coordsSize, cell);
+    coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 #endif
 
     // Dimensionalize coordinates for querying

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc	2013-05-08 22:30:43 UTC (rev 22004)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc	2013-05-09 00:00:30 UTC (rev 22005)
@@ -144,8 +144,6 @@
   PetscScalar* propertiesArray = propertiesVisitor.localArray();
 
   topology::CoordsVisitor coordsVisitor(dmMesh);
-  PetscScalar* coordsArray = NULL;
-  PetscInt coordsSize = 0;
 
   // Create arrays for querying.
   const int numDBProperties = _metadata.numDBProperties();
@@ -200,9 +198,11 @@
   for(PetscInt c = 0; c < numCells; ++c) {
     const PetscInt cell = cells[c];
     // Compute geometry information for current cell
-    coordsVisitor.getClosure(&coordsArray, &coordsSize, cell);
-    quadrature->computeGeometry(coordsArray, coordsSize, cell);
-    coordsVisitor.restoreClosure(&coordsArray, &coordsSize, cell);
+    PetscScalar* coordsCell = NULL;
+    PetscInt coordsSize = 0;
+    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
+    quadrature->computeGeometry(coordsCell, coordsSize, cell);
+    coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 
     const scalar_array& quadPtsNonDim = quadrature->quadPts();
     quadPtsGlobal = quadPtsNonDim;