[cig-commits] r21563 - in short/3D/PyLith/trunk: libsrc/pylith/materials unittests/libtests/materials unittests/libtests/materials/data

brad at geodynamics.org brad at geodynamics.org
Mon Mar 18 16:41:29 PDT 2013


Author: brad
Date: 2013-03-18 16:41:28 -0700 (Mon, 18 Mar 2013)
New Revision: 21563

Modified:
   short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.cc
Log:
Code cleanup.

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc	2013-03-18 23:41:14 UTC (rev 21562)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc	2013-03-18 23:41:28 UTC (rev 21563)
@@ -22,6 +22,9 @@
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/Stratum.hh" // USES StratumIS
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/utils/array.hh" // USES scalar_array, std::vector
 
@@ -93,7 +96,7 @@
 void
 pylith::materials::Material::normalizer(const spatialdata::units::Nondimensional& dim)
 { // normalizer
-  if (0 == _normalizer)
+  if (!_normalizer)
     _normalizer = new spatialdata::units::Nondimensional(dim);
   else
     *_normalizer = dim;
@@ -105,8 +108,8 @@
 pylith::materials::Material::initialize(const topology::Mesh& mesh,
 					feassemble::Quadrature<topology::Mesh>* quadrature)
 { // initialize
-  assert(0 != _dbProperties);
-  assert(0 != quadrature);
+  assert(_dbProperties);
+  assert(quadrature);
 
   // Get quadrature information
   const int numQuadPts = quadrature->numQuadPts();
@@ -114,49 +117,39 @@
   const int spaceDim = quadrature->spaceDim();
 
   // Get cells associated with material
-  DM              dmMesh = mesh.dmMesh();
-  IS              cellIS;
-  const PetscInt *cells;
-  PetscInt        numCells;
-  PetscErrorCode  err;
+  PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", _id);
+  const PetscInt numCells = materialIS.size();
+  const PetscInt* cells = materialIS.points();
 
-  assert(dmMesh);
-  err = DMPlexGetStratumIS(dmMesh, "material-id", _id, &cellIS);CHECK_PETSC_ERROR(err);
-  err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
-  assert(0 != cs);
+  const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();assert(cs);
 
   // Create field to hold physical properties.
-  delete _properties; _properties = new topology::Field<topology::Mesh>(mesh);
+  delete _properties; _properties = new topology::Field<topology::Mesh>(mesh);assert(_properties);
   _properties->label("properties");
-  assert(0 != _properties);
-  int fiberDim = numQuadPts * _numPropsQuadPt;
+  const int propsFiberDim = numQuadPts * _numPropsQuadPt;
   int_array cellsTmp(cells, numCells);
 
-  _properties->newSection(cellsTmp, fiberDim);
+  _properties->newSection(cellsTmp, propsFiberDim);
   _properties->allocate();
   _properties->zero();
-  PetscSection propertiesSection = _properties->petscSection();
-  PetscVec propertiesVec = _properties->localVector();
+  topology::VecVisitorMesh propertiesVisitor(*_properties);
+  PetscScalar* propertiesArray = propertiesVisitor.localArray();
 
 #if !defined(PRECOMPUTE_GEOMETRY)
-  scalar_array coordinatesCell(numBasis*spaceDim);
-  PetscSection coordSection;
-  PetscVec coordVec;
-  err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
-  err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
-  assert(coordSection);assert(coordVec);
+  topology::CoordsVisitor coordsVisitor(dmMesh);
+  PetscScalar* coordsArray = NULL;
+  PetscInt coordsSize = 0;
 #endif
 
   // Create arrays for querying.
   const int numDBProperties = _metadata.numDBProperties();
   scalar_array quadPtsGlobal(numQuadPts*spaceDim);
   scalar_array propertiesQuery(numDBProperties);
-  scalar_array propertiesCell(numQuadPts*_numPropsQuadPt);
+  scalar_array propertiesCell(propsFiberDim);
 
   // Setup database for quering for physical properties
-  assert(0 != _dbProperties);
+  assert(_dbProperties);
   _dbProperties->open();
   _dbProperties->queryVals(_metadata.dbProperties(),
 			   _metadata.numDBProperties());
@@ -164,32 +157,31 @@
   // Create field to hold state variables. We create the field even
   // if there is no initial state, because this we will use this field
   // to hold the state variables.
-  PetscSection stateVarsSection = PETSC_NULL;
-  PetscVec stateVarsVec = PETSC_NULL;
-  delete _stateVars; _stateVars = new topology::Field<topology::Mesh>(mesh);
+  delete _stateVars; _stateVars = new topology::Field<topology::Mesh>(mesh);assert(_stateVars);
   _stateVars->label("state variables");
-  fiberDim = numQuadPts * _numVarsQuadPt;
-  if (fiberDim > 0) {
-    assert(0 != _stateVars);
-    assert(0 != _properties);
-    _stateVars->newSection(*_properties, fiberDim);
+  const int stateVarsFiberDim = numQuadPts * _numVarsQuadPt;
+  topology::VecVisitorMesh* stateVarsVisitor = 0;
+  PetscScalar* stateVarsArray = NULL;
+  if (stateVarsFiberDim > 0) {
+    assert(_stateVars);
+    assert(_properties);
+    _stateVars->newSection(*_properties, stateVarsFiberDim);
     _stateVars->allocate();
     _stateVars->zero();
-    stateVarsSection = _stateVars->petscSection();
-    stateVarsVec     = _stateVars->localVector();
-    assert(stateVarsSection);assert(stateVarsVec);
+    stateVarsVisitor = new topology::VecVisitorMesh(*_stateVars);
+    stateVarsArray = stateVarsVisitor->localArray();
   } // if
 
+
   // Create arrays for querying
   const int numDBStateVars = _metadata.numDBStateVars();
   scalar_array stateVarsQuery;
   scalar_array stateVarsCell;
-  if (0 != _dbInitialState) {
-    assert(stateVarsSection);
+  if (_dbInitialState) {
     assert(numDBStateVars > 0);
     assert(_numVarsQuadPt > 0);
     stateVarsQuery.resize(numDBStateVars);
-    stateVarsCell.resize(numQuadPts*numDBStateVars);
+    stateVarsCell.resize(stateVarsFiberDim);
     
     // Setup database for querying for initial state variables
     _dbInitialState->open();
@@ -197,55 +189,41 @@
 			       _metadata.numDBStateVars());
   } // if
 
-  assert(0 != _normalizer);
+  assert(_normalizer);
   const PylithScalar lengthScale = _normalizer->lengthScale();
-  PetscScalar *propertiesArray, *stateVarsArray;
 
-  err = VecGetArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
-  if (stateVarsVec) {err = VecGetArray(stateVarsVec,  &stateVarsArray);CHECK_PETSC_ERROR(err);}
   for(PetscInt c = 0; c < numCells; ++c) {
     const PetscInt cell = cells[c];
     // Compute geometry information for current cell
 #if defined(PRECOMPUTE_GEOMETRY)
     quadrature->retrieveGeometry(cell);
 #else
-    PetscScalar *coords;
-    PetscInt     coordsSize;
-    err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
-    for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
-    quadrature->computeGeometry(coordinatesCell, cell);
-    err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+    coordsVisitor.getClosure(&coordsArray, &coordsSize, cell);
+    quadrature->computeGeometry(coordsArray, coordsSize, cell);
+    coordsVisitor.restoreClosure(&coordsArray, &coordsSize, cell);
 #endif
 
     const scalar_array& quadPtsNonDim = quadrature->quadPts();
     quadPtsGlobal = quadPtsNonDim;
-    _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
-				lengthScale);
+    _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(), lengthScale);
 
     // Loop over quadrature points in cell and query database
-    for (int iQuadPt=0, index=0; 
-	 iQuadPt < numQuadPts; 
-	 ++iQuadPt, index+=spaceDim) {
-      int err = _dbProperties->query(&propertiesQuery[0], numDBProperties,
-				     &quadPtsGlobal[index], spaceDim, cs);
+    for (int iQuadPt=0, index=0; iQuadPt < numQuadPts; ++iQuadPt, index+=spaceDim) {
+      int err = _dbProperties->query(&propertiesQuery[0], numDBProperties, &quadPtsGlobal[index], spaceDim, cs);
       if (err) {
 	std::ostringstream msg;
-	msg << "Could not find parameters for physical properties at \n"
-	    << "(";
+	msg << "Could not find parameters for physical properties at \n" << "(";
 	for (int i=0; i < spaceDim; ++i)
 	  msg << "  " << quadPtsGlobal[index+i];
 	msg << ") in material " << _label << "\n"
 	    << "using spatial database '" << _dbProperties->label() << "'.";
 	throw std::runtime_error(msg.str());
       } // if
-      _dbToProperties(&propertiesCell[iQuadPt*_numPropsQuadPt], 
-		      propertiesQuery);
-      _nondimProperties(&propertiesCell[iQuadPt*_numPropsQuadPt],
-			_numPropsQuadPt);
+      _dbToProperties(&propertiesCell[iQuadPt*_numPropsQuadPt], propertiesQuery);
+      _nondimProperties(&propertiesCell[iQuadPt*_numPropsQuadPt], _numPropsQuadPt);
 
-      if (0 != _dbInitialState) {
-	err = _dbInitialState->query(&stateVarsQuery[0], numDBStateVars,
-				     &quadPtsGlobal[index], spaceDim, cs);
+      if (_dbInitialState) {
+	err = _dbInitialState->query(&stateVarsQuery[0], numDBStateVars, &quadPtsGlobal[index], spaceDim, cs);
 	if (err) {
 	  std::ostringstream msg;
 	  msg << "Could not find initial state variables at \n" << "(";
@@ -256,39 +234,32 @@
 	      << "'.";
 	  throw std::runtime_error(msg.str());
 	} // if
-	_dbToStateVars(&stateVarsCell[iQuadPt*_numVarsQuadPt], 
-		       stateVarsQuery);
-	_nondimStateVars(&stateVarsCell[iQuadPt*_numVarsQuadPt],
-			 _numVarsQuadPt);
+	_dbToStateVars(&stateVarsCell[iQuadPt*_numVarsQuadPt], stateVarsQuery);
+	_nondimStateVars(&stateVarsCell[iQuadPt*_numVarsQuadPt], _numVarsQuadPt);
       } // if
 
     } // for
     // Insert cell contribution into fields
-    PetscInt dof, off, d;
-
-    err = PetscSectionGetDof(propertiesSection, cell, &dof);CHECK_PETSC_ERROR(err);
-    err = PetscSectionGetOffset(propertiesSection, cell, &off);CHECK_PETSC_ERROR(err);
-    assert(dof == numQuadPts*_numPropsQuadPt);
-    for(PetscInt d = 0; d < dof; ++d) {
+    const PetscInt off = propertiesVisitor.sectionOffset(cell);
+    assert(propsFiberDim == propertiesVisitor.sectionDof(cell));
+    for(PetscInt d = 0; d < propsFiberDim; ++d) {
       propertiesArray[off+d] = propertiesCell[d];
-    }
-    if (0 != _dbInitialState) {
-      err = PetscSectionGetDof(stateVarsSection, cell, &dof);CHECK_PETSC_ERROR(err);
-      err = PetscSectionGetOffset(stateVarsSection, cell, &off);CHECK_PETSC_ERROR(err);
-      assert(dof == numQuadPts*numDBStateVars);
-      for(PetscInt d = 0; d < dof; ++d) {
+    } // for
+    if (_dbInitialState) {
+      assert(stateVarsVisitor);
+      assert(stateVarsArray);
+      const PetscInt off = stateVarsVisitor->sectionOffset(cell);
+      assert(stateVarsFiberDim == stateVarsVisitor->sectionDof(cell));
+      for(PetscInt d = 0; d < stateVarsFiberDim; ++d) {
         stateVarsArray[off+d] = stateVarsCell[d];
-      }
-    }
+      } // for
+    } // if
   } // for
-  err = VecRestoreArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
-  if (stateVarsVec) {err = VecRestoreArray(stateVarsVec,  &stateVarsArray);CHECK_PETSC_ERROR(err);}
-  err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
+  delete stateVarsVisitor; stateVarsVisitor = 0;
 
   // Close databases
   _dbProperties->close();
-  if (0 != _dbInitialState)
+  if (_dbInitialState)
     _dbInitialState->close();
 } // initialize
 
@@ -336,13 +307,10 @@
 pylith::materials::Material::getField(topology::Field<topology::Mesh> *field,
 				      const char* name) const
 { // getField
-  // Logging of allocation is handled by getField() caller since it
-  // manages the memory for the field argument.
+  assert(field);
+  assert(_properties);
+  assert(_stateVars);
 
-  assert(0 != field);
-  assert(0 != _properties);
-  assert(0 != _stateVars);
-
   int propertyIndex = -1;
   int stateVarIndex = -1;
   _findField(&propertyIndex, &stateVarIndex, name);
@@ -354,17 +322,11 @@
   } // else
 
   // Get cell information
-  DM              dmMesh = field->mesh().dmMesh();
-  IS              cellIS;
-  const PetscInt *cells;
-  PetscInt        numCells;
-  PetscErrorCode  err;
+  PetscDM dmMesh = field->mesh().dmMesh();assert(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", _id);
+  const PetscInt numCells = materialIS.size();
+  const PetscInt* cells = materialIS.points();
 
-  assert(dmMesh);
-  err = DMPlexGetStratumIS(dmMesh, "material-id", _id, &cellIS);CHECK_PETSC_ERROR(err);
-  err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-
   topology::FieldBase::VectorFieldEnum fieldType = topology::FieldBase::OTHER;
       
   if (propertyIndex >= 0) { // If field is a property
@@ -373,17 +335,16 @@
     for (int i=0; i < propertyIndex; ++i)
       propOffset += _metadata.getProperty(i).fiberDim;
     const int fiberDim = _metadata.getProperty(propertyIndex).fiberDim;
+    topology::VecVisitorMesh propertiesVisitor(*_properties);
+    PetscScalar* propertiesArray = propertiesVisitor.localArray();
 
     // Get properties section
-    PetscSection propertiesSection = _properties->petscSection();
-    Vec          propertiesVec     = _properties->localVector();
-    assert(propertiesSection);
     PetscInt totalPropsFiberDimLocal = 0;
     PetscInt totalPropsFiberDim = 0;
-    if (numCells > 0) {err = PetscSectionGetDof(propertiesSection, cells[0], &totalPropsFiberDimLocal);CHECK_PETSC_ERROR(err);}
-    MPI_Allreduce((void *) &totalPropsFiberDimLocal, 
-                  (void *) &totalPropsFiberDim, 1, 
-                  MPIU_INT, MPI_MAX, field->mesh().comm());
+    if (numCells > 0) {
+      totalPropsFiberDimLocal = propertiesVisitor.sectionDof(cells[0]);
+    } // if
+    MPI_Allreduce((void *) &totalPropsFiberDimLocal, (void *) &totalPropsFiberDim, 1, MPIU_INT, MPI_MAX, field->mesh().comm());
     assert(totalPropsFiberDim > 0);
     const int numPropsQuadPt = _numPropsQuadPt;
     const int numQuadPts = totalPropsFiberDim / numPropsQuadPt;
@@ -391,10 +352,10 @@
     const int totalFiberDim = numQuadPts * fiberDim;
 
     // Allocate buffer for property field if necessary.
-    PetscSection fieldSection    = field->petscSection();
-    bool         useCurrentField = PETSC_FALSE;
-    PetscInt     pStart, pEnd;
-
+    bool useCurrentField = false;
+    PetscSection fieldSection = field->petscSection();
+    PetscInt pStart, pEnd;
+    PetscErrorCode err;
     err = PetscSectionGetChart(fieldSection, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
     if (pEnd < 0) {
       err = DMPlexGetHeightStratum(dmMesh, 0, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
@@ -403,10 +364,10 @@
       // check fiber dimension
       PetscInt totalFiberDimCurrentLocal = 0;
       PetscInt totalFiberDimCurrent = 0;
-      if (numCells > 0) {err = PetscSectionGetDof(fieldSection, cells[0], &totalFiberDimCurrentLocal);CHECK_PETSC_ERROR(err);}
-      MPI_Allreduce((void *) &totalFiberDimCurrentLocal, 
-                    (void *) &totalFiberDimCurrent, 1, 
-                    MPIU_INT, MPI_MAX, field->mesh().comm());
+      if (numCells > 0) {
+	err = PetscSectionGetDof(fieldSection, cells[0], &totalFiberDimCurrentLocal);CHECK_PETSC_ERROR(err);
+      } // if
+      MPI_Allreduce((void *) &totalFiberDimCurrentLocal, (void *) &totalFiberDimCurrent, 1, MPIU_INT, MPI_MAX, field->mesh().comm());
       assert(totalFiberDimCurrent > 0);
       useCurrentField = totalFiberDim == totalFiberDimCurrent;
     } // if
@@ -419,32 +380,26 @@
     field->label(name);
     field->scale(1.0);
     fieldType = _metadata.getProperty(propertyIndex).fieldType;
+    topology::VecVisitorMesh fieldVisitor(*field);
+    PetscScalar* fieldArray = fieldVisitor.localArray();
 
     // Buffer for property at cell's quadrature points
     scalar_array propertiesCell(numPropsQuadPt);
 
     // Loop over cells
-    Vec          fieldVec = field->localVector();
-    PetscScalar *fieldArray, *propertiesArray;
-
-    err = VecGetArray(fieldVec,      &fieldArray);CHECK_PETSC_ERROR(err);
-    err = VecGetArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
     for(PetscInt c = 0; c < numCells; ++c) {
       const PetscInt cell = cells[c];
-      PetscInt       off, poff;
-   
-      err = PetscSectionGetOffset(fieldSection,      cell, &off);CHECK_PETSC_ERROR(err);
-      err = PetscSectionGetOffset(propertiesSection, cell, &poff);CHECK_PETSC_ERROR(err);
+
+      const PetscInt poff = propertiesVisitor.sectionOffset(cell);
+      const PetscInt foff = fieldVisitor.sectionOffset(cell);
       for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
         for (int i=0; i < numPropsQuadPt; ++i)
           propertiesCell[i] = propertiesArray[iQuad*numPropsQuadPt + poff+i];
         _dimProperties(&propertiesCell[0], numPropsQuadPt);
         for (int i=0; i < fiberDim; ++i)
-          fieldArray[iQuad*fiberDim + off+i] = propertiesCell[propOffset+i];
+          fieldArray[iQuad*fiberDim + foff+i] = propertiesCell[propOffset+i];
       } // for
     } // for
-    err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
-    err = VecRestoreArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
   } else { // field is a state variable
     assert(stateVarIndex >= 0);
     
@@ -454,16 +409,16 @@
       varOffset += _metadata.getStateVar(i).fiberDim;
     const int fiberDim = _metadata.getStateVar(stateVarIndex).fiberDim;
 
-    // Get state variables section
-    PetscSection stateVarsSection = _stateVars->petscSection();
-    Vec          stateVarsVec     = _stateVars->localVector();
-    assert(stateVarsSection);assert(stateVarsVec);
+    // Get state variables
+    topology::VecVisitorMesh stateVarsVisitor(*_stateVars);
+    PetscScalar* stateVarsArray = stateVarsVisitor.localArray();
+
     PetscInt totalVarsFiberDimLocal = 0;
     PetscInt totalVarsFiberDim = 0;
-    if (numCells > 0) {err = PetscSectionGetDof(stateVarsSection, cells[0], &totalVarsFiberDimLocal);CHECK_PETSC_ERROR(err);}
-    MPI_Allreduce((void *) &totalVarsFiberDimLocal, 
-                  (void *) &totalVarsFiberDim, 1, 
-                  MPIU_INT, MPI_MAX, field->mesh().comm());
+    if (numCells > 0) {
+      totalVarsFiberDimLocal = stateVarsVisitor.sectionDof(cells[0]);
+    } // if
+    MPI_Allreduce((void*) &totalVarsFiberDimLocal, (void*) &totalVarsFiberDim, 1, MPIU_INT, MPI_MAX, field->mesh().comm());
     assert(totalVarsFiberDim > 0);
     const int numVarsQuadPt = _numVarsQuadPt;
     const int numQuadPts = totalVarsFiberDim / numVarsQuadPt;
@@ -471,17 +426,16 @@
     const int totalFiberDim = numQuadPts * fiberDim;
 
     // Allocate buffer for state variable field if necessary.
-    PetscSection fieldSection    = field->petscSection();
-    PetscVec          fieldVec        = field->localVector();
-    bool useCurrentField = fieldSection != PETSC_NULL;
+    PetscSection fieldSection = field->petscSection();
+    bool useCurrentField = fieldSection != NULL;
     if (fieldSection) {
       // check fiber dimension
       PetscInt totalFiberDimCurrentLocal = 0;
       PetscInt totalFiberDimCurrent = 0;
-      if (numCells > 0) {err = PetscSectionGetDof(fieldSection, cells[0], &totalFiberDimCurrentLocal);CHECK_PETSC_ERROR(err);}
-      MPI_Allreduce((void *) &totalFiberDimCurrentLocal, 
-                    (void *) &totalFiberDimCurrent, 1, 
-                    MPIU_INT, MPI_MAX, field->mesh().comm());
+      if (numCells > 0) {
+	PetscErrorCode err = PetscSectionGetDof(fieldSection, cells[0], &totalFiberDimCurrentLocal);CHECK_PETSC_ERROR(err);
+      } // if
+      MPI_Allreduce((void*) &totalFiberDimCurrentLocal, (void*) &totalFiberDimCurrent, 1, MPIU_INT, MPI_MAX, field->mesh().comm());
       assert(totalFiberDimCurrent > 0);
       useCurrentField = totalFiberDim == totalFiberDimCurrent;
     } // if
@@ -494,34 +448,28 @@
     fieldType = _metadata.getStateVar(stateVarIndex).fieldType;
     field->label(name);
     field->scale(1.0);
+    topology::VecVisitorMesh fieldVisitor(*field);
+    PetscScalar* fieldArray = fieldVisitor.localArray();
 
     // Buffer for state variable at cell's quadrature points
     scalar_array stateVarsCell(numVarsQuadPt);
     
     // Loop over cells
-    PetscScalar  *fieldArray, *stateVarsArray;
-
-    err = VecGetArray(fieldVec,     &fieldArray);CHECK_PETSC_ERROR(err);
-    err = VecGetArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);
     for(PetscInt c = 0; c < numCells; ++c) {
       const PetscInt cell = cells[c];
-      PetscInt       off, soff;
-      
-      err = PetscSectionGetOffset(fieldSection,     cell, &off);CHECK_PETSC_ERROR(err);
-      err = PetscSectionGetOffset(stateVarsSection, cell, &soff);CHECK_PETSC_ERROR(err);
+
+      const PetscInt foff = fieldVisitor.sectionOffset(cell);
+      const PetscInt soff = stateVarsVisitor.sectionOffset(cell);
       for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
-        for (int i=0; i < numVarsQuadPt; ++i)
+        for (int i=0; i < numVarsQuadPt; ++i) {
           stateVarsCell[i] = stateVarsArray[iQuad*numVarsQuadPt + soff+i];
-        _dimStateVars(&stateVarsCell[0], numVarsQuadPt);
+	} // for
+	_dimStateVars(&stateVarsCell[0], numVarsQuadPt);
         for (int i=0; i < fiberDim; ++i)
-          fieldArray[iQuad*fiberDim + off+i] = stateVarsCell[varOffset+i];
+          fieldArray[iQuad*fiberDim + foff+i] = stateVarsCell[varOffset+i];
       } // for
     } // for
-    err = VecRestoreArray(fieldVec,     &fieldArray);CHECK_PETSC_ERROR(err);
-    err = VecRestoreArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);
   } // if/else
-  err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
 
   topology::FieldBase::VectorFieldEnum multiType = 
     topology::FieldBase::MULTI_OTHER;
@@ -558,8 +506,8 @@
 					int* stateVarIndex,
 					const char* name) const
 { // _findField
-  assert(0 != propertyIndex);
-  assert(0 != stateVarIndex);
+  assert(propertyIndex);
+  assert(stateVarIndex);
 
   *propertyIndex = -1;
   *stateVarIndex = -1;

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc	2013-03-18 23:41:14 UTC (rev 21562)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc	2013-03-18 23:41:28 UTC (rev 21563)
@@ -25,6 +25,8 @@
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/Stratum.hh" // USES StratumIS
+#include "pylith/topology/VisitorMesh.hh" // USES VisitorMesh
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 #include "pylith/materials/ElasticIsotropic3D.hh" // USES ElasticIsotropic3D
 #include "pylith/materials/ElasticStrain1D.hh" // USES ElasticStrain1D
@@ -46,10 +48,6 @@
 CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestElasticMaterial );
 
 // ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
 // Test dbInitialStress()
 void
 pylith::materials::TestElasticMaterial::testDBInitialStress(void)
@@ -93,16 +91,10 @@
 
   // Get cells associated with material
   const int materialId = 24;
-  PetscDM dmMesh     = mesh.dmMesh();
-  PetscIS cellIS;
-  const PetscInt *cells;
-  PetscInt numCells;
-  PetscErrorCode err = 0;
-
-  assert(dmMesh);
-  err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
-  err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", materialId);
+  const PetscInt* cells = materialIS.points();
+  const PetscInt numCells = materialIS.size();
   PetscInt cell = cells[0];
 
   const PylithScalar tolerance = 1.0e-06;
@@ -110,39 +102,34 @@
   const int numQuadPts = data.numLocs;
 
   // Test initialStress field
-  CPPUNIT_ASSERT(material._initialFields);
-  PetscSection stressSection = material._initialFields->get("initial stress").petscSection();
-  Vec          stressVec     = material._initialFields->get("initial stress").localVector();
-  CPPUNIT_ASSERT(stressSection);CPPUNIT_ASSERT(stressVec);
-  PetscInt dof, off, fiberDim = numQuadPts * tensorSize;
-  err = PetscSectionGetDof(stressSection, cell, &dof);CHECK_PETSC_ERROR(err);
-  err = PetscSectionGetOffset(stressSection, cell, &off);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
   const PylithScalar *initialStressE = data.initialStress;
-  PetscScalar        *initialStress;
   CPPUNIT_ASSERT(initialStressE);
-  err = VecGetArray(stressVec, &initialStress);CHECK_PETSC_ERROR(err);
+  int fiberDim = numQuadPts * tensorSize;
+
+  CPPUNIT_ASSERT(material._initialFields);
+  topology::Field<topology::Mesh>& stressField = material._initialFields->get("initial stress");
+  topology::VecVisitorMesh stressVisitor(stressField);
+  const PetscScalar* initialStress = stressVisitor.localArray();
+  PetscInt off = stressVisitor.sectionOffset(cell);
+  CPPUNIT_ASSERT_EQUAL(fiberDim, stressVisitor.sectionDof(cell));
   for(int i = 0; i < fiberDim; ++i) {
     CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStress[off+i]/initialStressE[i]*data.pressureScale, tolerance);
-  }
-  err = VecRestoreArray(stressVec, &initialStress);CHECK_PETSC_ERROR(err);
+  } // for
 
   // Test initialStrain field
-  CPPUNIT_ASSERT(material._initialFields);
-  PetscSection strainSection = material._initialFields->get("initial strain").petscSection();
-  Vec          strainVec     = material._initialFields->get("initial strain").localVector();
-  fiberDim = numQuadPts * tensorSize;
-  err = PetscSectionGetDof(strainSection, cell, &dof);CHECK_PETSC_ERROR(err);
-  err = PetscSectionGetOffset(strainSection, cell, &off);CHECK_PETSC_ERROR(err);
-  CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
   const PylithScalar *initialStrainE = data.initialStrain;
-  PetscScalar        *initialStrain;
   CPPUNIT_ASSERT(initialStrainE);
-  err = VecGetArray(strainVec, &initialStrain);CHECK_PETSC_ERROR(err);
+  fiberDim = numQuadPts * tensorSize;
+
+  CPPUNIT_ASSERT(material._initialFields);
+  topology::Field<topology::Mesh>& strainField = material._initialFields->get("initial strain");
+  topology::VecVisitorMesh strainVisitor(strainField);
+  const PetscScalar* initialStrain = strainVisitor.localArray();
+  off = strainVisitor.sectionOffset(cell);
+  CPPUNIT_ASSERT_EQUAL(fiberDim, strainVisitor.sectionDof(cell));
   for(int i = 0; i < fiberDim; ++i) {
     CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStrain[off+i]/initialStrainE[i], tolerance);
-  }
-  err = VecRestoreArray(strainVec, &initialStrain);CHECK_PETSC_ERROR(err);
+  } // for
 
   // Test cell arrays
   size_t size = data.numLocs*data.numPropsQuadPt;
@@ -180,8 +167,6 @@
     } // switch
   size = data.numLocs*numElasticConsts;
   CPPUNIT_ASSERT_EQUAL(size, material._elasticConstsCell.size());
-  err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
 } // testInitialize
 
 // ----------------------------------------------------------------------
@@ -196,19 +181,11 @@
 
   // Get cells associated with material
   const int materialId = 24;
-  PetscDM dmMesh     = mesh.dmMesh();
-  PetscIS cellIS;
-  const PetscInt *cells;
-  PetscInt numCells;
-  PetscErrorCode err = 0;
-
-  assert(dmMesh);
-  err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
-  err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", materialId);
+  const PetscInt* cells = materialIS.points();
+  const PetscInt numCells = materialIS.size();
   PetscInt cell = cells[0];
-  err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
 
   material.retrievePropsAndVars(cell);
 
@@ -268,19 +245,11 @@
 
   // Get cells associated with material
   const int materialId = 24;
-  PetscDM dmMesh = mesh.dmMesh();
-  PetscIS cellIS;
-  const PetscInt *cells;
-  PetscInt numCells;
-  PetscErrorCode err = 0;
-
-  assert(dmMesh);
-  err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
-  err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", materialId);
+  const PetscInt* cells = materialIS.points();
+  const PetscInt numCells = materialIS.size();
   PetscInt cell = cells[0];
-  err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
 
   material.retrievePropsAndVars(cell);
   const scalar_array& density = material.calcDensity();
@@ -309,19 +278,11 @@
 
   // Get cells associated with material
   const int materialId = 24;
-  PetscDM dmMesh = mesh.dmMesh();
-  PetscIS cellIS;
-  const PetscInt *cells;
-  PetscInt numCells;
-  PetscErrorCode err = 0;
-
-  assert(dmMesh);
-  err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
-  err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", materialId);
+  const PetscInt* cells = materialIS.points();
+  const PetscInt numCells = materialIS.size();
   PetscInt cell = cells[0];
-  err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
 
   const int tensorSize = material._tensorSize;
   const int numQuadPts = data.numLocs;
@@ -353,19 +314,11 @@
 
   // Get cells associated with material
   const int materialId = 24;
-  PetscDM dmMesh = mesh.dmMesh();
-  PetscIS cellIS;
-  const PetscInt *cells;
-  PetscInt numCells;
-  PetscErrorCode  err;
-
-  assert(dmMesh);
-  err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
-  err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", materialId);
+  const PetscInt* cells = materialIS.points();
+  const PetscInt numCells = materialIS.size();
   PetscInt cell = cells[0];
-  err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
 
   const int tensorSize = material._tensorSize;
   const int numQuadPts = data.numLocs;
@@ -423,19 +376,11 @@
 
   // Get cells associated with material
   const int materialId = 24;
-  PetscDM dmMesh     = mesh.dmMesh();
-  PetscIS cellIS;
-  const PetscInt *cells;
-  PetscInt numCells;
-  PetscErrorCode err = 0;
-
-  assert(dmMesh);
-  err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
-  err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", materialId);
+  const PetscInt* cells = materialIS.points();
+  const PetscInt numCells = materialIS.size();
   PetscInt cell = cells[0];
-  err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
 
   material.retrievePropsAndVars(cell);
   const PylithScalar dt = material.stableTimeStepImplicit(mesh);
@@ -457,11 +402,11 @@
 
   // Get cells associated with material
   const int materialId = 24;
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  assert(!sieveMesh.isNull());
-  const ALE::Obj<SieveMesh::label_sequence>& cells = 
-    sieveMesh->getLabelStratum("material-id", materialId);
-  SieveMesh::point_type cell = *cells->begin();
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", materialId);
+  const PetscInt* cells = materialIS.points();
+  const PetscInt numCells = materialIS.size();
+  PetscInt cell = cells[0];
 
   // Setup quadrature
   feassemble::Quadrature<topology::Mesh> quadrature;
@@ -794,10 +739,9 @@
 // ----------------------------------------------------------------------
 // Setup mesh and material.
 void
-pylith::materials::TestElasticMaterial::_initialize(
-					  topology::Mesh* mesh,
-					  ElasticStrain1D* material,
-					  const ElasticStrain1DData* data)
+pylith::materials::TestElasticMaterial::_initialize(topology::Mesh* mesh,
+						    ElasticStrain1D* material,
+						    const ElasticStrain1DData* data)
 { // _initialize
   CPPUNIT_ASSERT(mesh);
   CPPUNIT_ASSERT(material);

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc	2013-03-18 23:41:14 UTC (rev 21562)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc	2013-03-18 23:41:28 UTC (rev 21563)
@@ -26,6 +26,8 @@
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/Stratum.hh" // USES StratumIS
+#include "pylith/topology/VisitorMesh.hh" // USES VisitorMesh
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 #include "pylith/materials/ElasticIsotropic3D.hh" // USES ElasticIsotropic3D
 #include "pylith/materials/ElasticStrain1D.hh" // USES ElasticStrain1D
@@ -47,10 +49,6 @@
 CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestMaterial );
 
 // ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
 // Test id()
 void
 pylith::materials::TestMaterial::testID(void)
@@ -98,7 +96,7 @@
   ElasticIsotropic3D material;
   material.dbProperties(&db);
   
-  CPPUNIT_ASSERT(0 != material._dbProperties);
+  CPPUNIT_ASSERT(material._dbProperties);
   CPPUNIT_ASSERT_EQUAL(label, std::string(material._dbProperties->label()));
 } // testDBProperties
 
@@ -114,7 +112,7 @@
   ElasticIsotropic3D material;
   material.dbInitialState(&db);
   
-  CPPUNIT_ASSERT(0 != material._dbInitialState);
+  CPPUNIT_ASSERT(material._dbInitialState);
   CPPUNIT_ASSERT_EQUAL(label, std::string(material._dbInitialState->label()));
 } // testDBStateVars
 
@@ -130,7 +128,7 @@
   ElasticIsotropic3D material;
   material.normalizer(normalizer);
   
-  CPPUNIT_ASSERT(0 != material._normalizer);
+  CPPUNIT_ASSERT(material._normalizer);
   CPPUNIT_ASSERT_EQUAL(lengthScale, material._normalizer->lengthScale());
 } // testNormalizer
 
@@ -181,6 +179,18 @@
   cs.initialize();
   mesh.coordsys(&cs);
 
+  spatialdata::units::Nondimensional normalizer;
+  const PylithScalar lengthScale = 1.0e+3;
+  const PylithScalar pressureScale = 2.25e+10;
+  const PylithScalar timeScale = 2.0;
+  const PylithScalar velocityScale = lengthScale / timeScale;
+  const PylithScalar densityScale = pressureScale / (velocityScale*velocityScale);
+  normalizer.lengthScale(lengthScale);
+  normalizer.pressureScale(pressureScale);
+  normalizer.timeScale(timeScale);
+  normalizer.densityScale(densityScale);
+  mesh.nondimensionalize(normalizer);
+
   // Setup quadrature
   feassemble::Quadrature<topology::Mesh> quadrature;
   feassemble::GeometryLine1D geometry;
@@ -209,17 +219,11 @@
 
   // Get cells associated with material
   const PetscInt  materialId = 24;
-  DM              dmMesh     = mesh.dmMesh();
-  IS              cellIS;
-  const PetscInt *cells;
-  PetscInt        numCells;
-  PetscErrorCode  err;
+  PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+  topology::StratumIS materialIS(dmMesh, "material-id", materialId);
+  const PetscInt* cells = materialIS.points();
+  const PetscInt numCells = materialIS.size();
 
-  assert(dmMesh);
-  err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
-  err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
-  err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-
   // Compute geometry for cells
   quadrature.initializeGeometry();
 #if defined(PRECOMPUTE_GEOMETRY)
@@ -233,8 +237,6 @@
   db.ioHandler(&dbIO);
   db.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
   
-  spatialdata::units::Nondimensional normalizer;
-
   ElasticStrain1D material;
   material.dbProperties(&db);
   material.id(materialId);
@@ -256,44 +258,37 @@
   const PylithScalar muE[] = { muA, muB };
   const PylithScalar lambdaE[] = { lambdaA, lambdaB };
 
-  SieveMesh::point_type cell = cells[0];
+  PetscInt cell = cells[0];
   const PylithScalar tolerance = 1.0e-06;
 
-  CPPUNIT_ASSERT(0 != material._properties);
-  PetscSection propertiesSection = material._properties->petscSection();
-  Vec          propertiesVec     = material._properties->localVector();
-  CPPUNIT_ASSERT(propertiesSection);CPPUNIT_ASSERT(propertiesVec);
+  CPPUNIT_ASSERT(material._properties);
+  topology::VecVisitorMesh propertiesVisitor(*material._properties);
+  const PetscScalar* propertiesArray = propertiesVisitor.localArray();
+  const PetscInt off = propertiesVisitor.sectionOffset(cell);
 
   const int p_density = 0;
   const int p_mu = 1;
   const int p_lambda = 2;
-  PetscScalar *propertiesArray;
-  PetscInt     off;
 
-  err = VecGetArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
-  err = PetscSectionGetOffset(propertiesSection, cell, &off);CHECK_PETSC_ERROR(err);
   // density
   for (int i=0; i < numQuadPts; ++i) {
     const int index = i*material._numPropsQuadPt + p_density;
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesArray[off+index]/densityE[i],
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesArray[off+index]/densityE[i]*densityScale,
 				 tolerance);
   } // for
   
   // mu
   for (int i=0; i < numQuadPts; ++i) {
     const int index = i*material._numPropsQuadPt + p_mu;
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesArray[off+index]/muE[i], tolerance);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesArray[off+index]/muE[i]*pressureScale, tolerance);
   } // for
   
   // lambda
   for (int i=0; i < numQuadPts; ++i) {
     const int index = i*material._numPropsQuadPt + p_lambda;
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesArray[off+index]/lambdaE[i], 
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesArray[off+index]/lambdaE[i]*pressureScale, 
 				 tolerance);
   } // for
-  err = VecRestoreArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
-  err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-  err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
 } // testInitialize
 
 // ----------------------------------------------------------------------
@@ -353,8 +348,8 @@
 void
 pylith::materials::TestMaterial::testDBToProperties(void)
 { // testDBToProperties
-  CPPUNIT_ASSERT(0 != _material);
-  CPPUNIT_ASSERT(0 != _data);
+  CPPUNIT_ASSERT(_material);
+  CPPUNIT_ASSERT(_data);
   
   // Check to make sure names of Metadata values match names of test
   // data values (consistency check).
@@ -398,8 +393,8 @@
 void
 pylith::materials::TestMaterial::testNonDimProperties(void)
 { // testNonDimProperties
-  CPPUNIT_ASSERT(0 != _material);
-  CPPUNIT_ASSERT(0 != _data);
+  CPPUNIT_ASSERT(_material);
+  CPPUNIT_ASSERT(_data);
   
   const int numLocs = _data->numLocs;
   const int propertiesSize = _data->numPropsQuadPt;
@@ -414,7 +409,7 @@
     
     const PylithScalar* const propertiesNondimE =
       &_data->propertiesNondim[iLoc*propertiesSize];
-    CPPUNIT_ASSERT(0 != propertiesNondimE);
+    CPPUNIT_ASSERT(propertiesNondimE);
 
     const PylithScalar tolerance = 1.0e-06;
     for (int i=0; i < propertiesSize; ++i) {
@@ -434,8 +429,8 @@
 void
 pylith::materials::TestMaterial::testDimProperties(void)
 { // testDimProperties
-  CPPUNIT_ASSERT(0 != _material);
-  CPPUNIT_ASSERT(0 != _data);
+  CPPUNIT_ASSERT(_material);
+  CPPUNIT_ASSERT(_data);
   
   const int numLocs = _data->numLocs;
   const int propertiesSize = _data->numPropsQuadPt;
@@ -449,7 +444,7 @@
     
     const PylithScalar* const propertiesE =
       &_data->properties[iLoc*propertiesSize];
-    CPPUNIT_ASSERT(0 != propertiesE);
+    CPPUNIT_ASSERT(propertiesE);
 
     const PylithScalar tolerance = 1.0e-06;
     for (int i=0; i < propertiesSize; ++i) {
@@ -469,8 +464,8 @@
 void
 pylith::materials::TestMaterial::testDBToStateVars(void)
 { // testDBToStateVars
-  CPPUNIT_ASSERT(0 != _material);
-  CPPUNIT_ASSERT(0 != _data);
+  CPPUNIT_ASSERT(_material);
+  CPPUNIT_ASSERT(_data);
   
   // Check to make sure names of Metadata values match names of test
   // data values (consistency check).
@@ -497,7 +492,7 @@
     
     const PylithScalar* const stateVarsE = 
       (stateVarsSize > 0) ? &_data->stateVars[iLoc*stateVarsSize] : 0;
-    CPPUNIT_ASSERT( (0 < stateVarsSize && 0 != stateVarsE) ||
+    CPPUNIT_ASSERT( (0 < stateVarsSize && stateVarsE) ||
 		    (0 == stateVarsSize && 0 == stateVarsE) );
     const PylithScalar tolerance = 1.0e-06;
     for (int i=0; i < stateVarsSize; ++i) {
@@ -517,8 +512,8 @@
 void
 pylith::materials::TestMaterial::testNonDimStateVars(void)
 { // testNonDimStateVars
-  CPPUNIT_ASSERT(0 != _material);
-  CPPUNIT_ASSERT(0 != _data);
+  CPPUNIT_ASSERT(_material);
+  CPPUNIT_ASSERT(_data);
   
   const int numLocs = _data->numLocs;
   const int stateVarsSize = _data->numVarsQuadPt;
@@ -532,7 +527,7 @@
     
     const PylithScalar* const stateVarsNondimE =
       (stateVarsSize > 0) ? &_data->stateVarsNondim[iLoc*stateVarsSize] : 0;
-    CPPUNIT_ASSERT( (0 < stateVarsSize && 0 != stateVarsNondimE) ||
+    CPPUNIT_ASSERT( (0 < stateVarsSize && stateVarsNondimE) ||
 		    (0 == stateVarsSize && 0 == stateVarsNondimE) );
 
     const PylithScalar tolerance = 1.0e-06;
@@ -553,8 +548,8 @@
 void
 pylith::materials::TestMaterial::testDimStateVars(void)
 { // testDimStateVars
-  CPPUNIT_ASSERT(0 != _material);
-  CPPUNIT_ASSERT(0 != _data);
+  CPPUNIT_ASSERT(_material);
+  CPPUNIT_ASSERT(_data);
   
   const int numLocs = _data->numLocs;
   const int stateVarsSize = _data->numVarsQuadPt;
@@ -568,7 +563,7 @@
     
     const PylithScalar* const stateVarsE =
       (stateVarsSize > 0) ? &_data->stateVars[iLoc*stateVarsSize] : 0;
-    CPPUNIT_ASSERT( (0 < stateVarsSize && 0 != stateVarsE) ||
+    CPPUNIT_ASSERT( (0 < stateVarsSize && stateVarsE) ||
 		    (0 == stateVarsSize && 0 == stateVarsE) );
 
     const PylithScalar tolerance = 1.0e-06;

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.cc	2013-03-18 23:41:14 UTC (rev 21562)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.cc	2013-03-18 23:41:28 UTC (rev 21563)
@@ -39,11 +39,13 @@
   stateVars(0),
   propertiesNondim(0),
   stateVarsNondim(0),
-  lengthScale(0),
-  timeScale(0),
-  pressureScale(0),
+  lengthScale(1.0e+3),
+  timeScale(2.0),
+  pressureScale(2.25e+10),
   densityScale(0)
 { // constructor
+  const PylithScalar velScale = lengthScale / timeScale;
+  densityScale = pressureScale / (velScale*velScale);
 } // constructor
 
 // ----------------------------------------------------------------------



More information about the CIG-COMMITS mailing list