[cig-commits] r6476 - in short/3D/PyLith/trunk: libsrc/feassemble libsrc/materials libsrc/meshio modulesrc/feassemble modulesrc/meshio modulesrc/topology unittests/libtests/feassemble unittests/libtests/materials unittests/libtests/meshio

knepley at geodynamics.org knepley at geodynamics.org
Fri Mar 30 16:49:17 PDT 2007


Author: knepley
Date: 2007-03-30 16:49:16 -0700 (Fri, 30 Mar 2007)
New Revision: 6476

Removed:
   short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc
   short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh
Modified:
   short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
   short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
   short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh
   short/3D/PyLith/trunk/libsrc/feassemble/ParameterManager.cc
   short/3D/PyLith/trunk/libsrc/feassemble/ParameterManager.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.hh
   short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
   short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
   short/3D/PyLith/trunk/libsrc/materials/Material.cc
   short/3D/PyLith/trunk/libsrc/materials/Material.hh
   short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
   short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh
   short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
   short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src
   short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
   short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.hh
Log:
Update to latest Sieve


Modified: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -61,6 +61,7 @@
 // finite elements.
 void
 pylith::feassemble::ExplicitElasticity::integrateConstant(
+                  const ALE::Obj<Mesh>& m,
 			      const ALE::Obj<real_section_type>& b,
 			      const ALE::Obj<real_section_type>& dispT,
 			      const ALE::Obj<real_section_type>& dispTmdt,
@@ -69,11 +70,8 @@
   assert(0 != _quadrature);
 
   // Get information about section
-  const topology_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = dispT->getTopology();
-  const ALE::Obj<topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
-  const topology_type::label_sequence::iterator cellsEnd = cells->end();
+  const ALE::Obj<Mesh::label_sequence>& cells    = m->heightStratum(0);
+  const Mesh::label_sequence::iterator  cellsEnd = cells->end();
 
   // Get parameters used in integration.
   const double dt = _dt;
@@ -89,20 +87,20 @@
   const int spaceDim = _quadrature->spaceDim();
   const int cellDim = _quadrature->cellDim();
 
-  for (topology_type::label_sequence::iterator cellIter=cells->begin();
+  for (Mesh::label_sequence::iterator cellIter=cells->begin();
        cellIter != cellsEnd;
        ++cellIter) {
     // Compute geometry information for current cell
-    _quadrature->computeGeometry(coordinates, *cellIter);
+    _quadrature->computeGeometry(m, coordinates, *cellIter);
 
     // Reset element vector to zero
     _resetCellVector();
 
     // Restrict input fields to cell
     const real_section_type::value_type* dispTCell = 
-      dispT->restrict(patch, *cellIter);
+      m->restrict(dispT, *cellIter);
     const real_section_type::value_type* dispTmdtCell = 
-      dispTmdt->restrict(patch, *cellIter);
+      m->restrict(dispTmdt, *cellIter);
 
     // Get cell geometry information that depends on cell
     const double* basis = _quadrature->basis();
@@ -110,7 +108,7 @@
     const double* jacobianDet = _quadrature->jacobianDet();
 
     // Get material physical properties at quadrature points for this cell
-    _material->calcProperties(*cellIter, patch, numQuadPts);
+    _material->calcProperties(*cellIter, numQuadPts);
     const double* density = _material->density();
     const double* elasticConsts = _material->elasticConsts();
     const int numElasticConsts = _material->numElasticConsts();
@@ -291,7 +289,7 @@
     } // if/else
 
     // Assemble cell contribution into field
-    b->updateAdd(patch, *cellIter, _cellVector);
+    m->updateAdd(b, *cellIter, _cellVector);
   } // for
 } // integrateConstant
 
@@ -300,17 +298,15 @@
 void
 pylith::feassemble::ExplicitElasticity::integrateJacobian(
 			     PetscMat* mat,
+                 const ALE::Obj<Mesh>& m,
 			     const ALE::Obj<real_section_type>& dispT,
 			     const ALE::Obj<real_section_type>& coordinates)
 { // integrateJacobian
   assert(0 != _quadrature);
 
   // Get information about section
-  const topology_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = dispT->getTopology();
-  const ALE::Obj<topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
-  const topology_type::label_sequence::iterator cellsEnd = cells->end();
+  const ALE::Obj<Mesh::label_sequence>& cells    = m->heightStratum(0);
+  const Mesh::label_sequence::iterator  cellsEnd = cells->end();
 
   // Get parameters used in integration.
   const double dt = _dt;
@@ -324,11 +320,11 @@
   const double* quadWts = _quadrature->quadWts();
   const int numBasis = _quadrature->numCorners();
   const int spaceDim = _quadrature->spaceDim();
-  for (topology_type::label_sequence::iterator cellIter=cells->begin();
+  for (Mesh::label_sequence::iterator cellIter=cells->begin();
        cellIter != cellsEnd;
        ++cellIter) {
     // Compute geometry information for current cell
-    _quadrature->computeGeometry(coordinates, *cellIter);
+    _quadrature->computeGeometry(m, coordinates, *cellIter);
 
     // Reset element vector to zero
     _resetCellMatrix();
@@ -338,7 +334,7 @@
     const double* jacobianDet = _quadrature->jacobianDet();
 
     // Get material physical properties at quadrature points for this cell
-    _material->calcProperties(*cellIter, patch, numQuadPts);
+    _material->calcProperties(*cellIter, numQuadPts);
     const double* density = _material->density();
 
     // Compute Jacobian for cell
@@ -364,7 +360,9 @@
       throw std::runtime_error("Logging PETSc flops failed.");
     
     // Assemble cell contribution into field
-    err = assembleMatrix(*mat, *cellIter, _cellMatrix, ADD_VALUES);
+    const ALE::Obj<Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, "default", dispT->getAtlas());
+
+    err = updateOperator(*mat, m, dispT, globalOrder, *cellIter, _cellMatrix, ADD_VALUES);
     if (err)
       throw std::runtime_error("Update to PETSc Mat failed.");
   } // for

Modified: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -106,7 +106,8 @@
    * @param dispTmdt Displacement field at time t-dt
    * @param coordinates Field of cell vertex coordinates
    */
-  void integrateConstant(const ALE::Obj<real_section_type>& b,
+  void integrateConstant(const ALE::Obj<Mesh>& m,
+             const ALE::Obj<real_section_type>& b,
 			 const ALE::Obj<real_section_type>& dispT,
 			 const ALE::Obj<real_section_type>& dispTmdt,
 			 const ALE::Obj<real_section_type>& coordinates);
@@ -118,6 +119,7 @@
    * @param coordinates Field of cell vertex coordinates
    */
   void integrateJacobian(PetscMat* mat,
+             const ALE::Obj<Mesh>& m,
 			 const ALE::Obj<real_section_type>& dispT,
 			 const ALE::Obj<real_section_type>& coordinates);
   

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -43,9 +43,8 @@
 // PUBLIC TYPEDEFS //////////////////////////////////////////////////////
 public :
 
-  typedef ALE::Mesh Mesh;
-  typedef Mesh::topology_type topology_type;
-  typedef topology_type::point_type point_type;
+  typedef ALE::Field::Mesh        Mesh;
+  typedef ALE::Mesh::point_type   point_type;
   typedef Mesh::real_section_type real_section_type;
 
 // PUBLIC MEMBERS ///////////////////////////////////////////////////////

Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -82,7 +82,8 @@
    * @param coordinates Field of cell vertex coordinates
    */
   virtual 
-  void integrateConstant(const ALE::Obj<real_section_type>& fieldOut,
+  void integrateConstant(const ALE::Obj<Mesh>& mesh,
+             const ALE::Obj<real_section_type>& fieldOut,
 			 const ALE::Obj<real_section_type>& fieldInT,
 			 const ALE::Obj<real_section_type>& fieldInTmdt,
 			 const ALE::Obj<real_section_type>& coordinates) = 0;
@@ -95,6 +96,7 @@
    */
   virtual 
   void integrateJacobian(PetscMat* mat,
+             const ALE::Obj<Mesh>& mesh,
 			 const ALE::Obj<real_section_type>& fieldIn,
 			 const ALE::Obj<real_section_type>& coordinates) = 0;
 

Deleted: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -1,330 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-//                           Brad T. Aagaard
-//                        U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#include <portinfo>
-
-#include "IntegratorInertia.hh" // implementation of class methods
-
-#include "Quadrature.hh" // USES Quadrature
-
-#include "petscmat.h" // USES PetscMat
-#include "spatialdata/spatialdb/SpatialDB.hh"
-
-#include <assert.h> // USES assert()
-#include <stdexcept> // USES std::runtime_error
-
-// ----------------------------------------------------------------------
-// Constructor
-pylith::feassemble::IntegratorInertia::IntegratorInertia(void) :
-  Integrator()
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor
-pylith::feassemble::IntegratorInertia::~IntegratorInertia(void)
-{ // destructor
-} // destructor
-  
-// ----------------------------------------------------------------------
-// Copy constructor.
-pylith::feassemble::IntegratorInertia::IntegratorInertia(const IntegratorInertia& i) :
-  Integrator(i)
-{ // copy constructor
-} // copy constructor
-
-// ----------------------------------------------------------------------
-// Integrate inertial term for 3-D finite elements.
-void
-pylith::feassemble::IntegratorInertia::integrateAction(
-			      const ALE::Obj<real_section_type>& fieldOut,
-			      const ALE::Obj<real_section_type>& fieldIn,
-			      const ALE::Obj<real_section_type>& coordinates)
-{ // integrateAction
-  assert(0 != _quadrature);
-  assert(0 != _density.isNull());
-
-  // Get information about section
-  const topology_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = fieldIn->getTopology();
-  const ALE::Obj<topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
-  const topology_type::label_sequence::iterator cellsEnd = cells->end();
-
-  // Allocate vector for cell values (if necessary)
-  _initCellVector();
-
-  for (topology_type::label_sequence::iterator cellIter=cells->begin();
-       cellIter != cellsEnd;
-       ++cellIter) {
-    // Compute geometry information for current cell
-    _quadrature->computeGeometry(coordinates, *cellIter);
-
-    // Reset element vector to zero
-    _resetCellVector();
-
-    // Restrict input field to cell
-    const real_section_type::value_type* fieldInCell = 
-      fieldIn->restrict(patch, *cellIter);
-
-    // Get cell geometry information
-    const int numQuadPts = _quadrature->numQuadPts();
-    const double* basis = _quadrature->basis();
-    const double* quadPts = _quadrature->quadPts();
-    const double* quadWts = _quadrature->quadWts();
-    const double* jacobianDet = _quadrature->jacobianDet();
-    const int numBasis = _quadrature->numCorners();
-    const int spaceDim = _quadrature->spaceDim();
-
-    // Restrict density from material database to quadrature points for this cell
-    const real_section_type::value_type* density = 
-      _density->restrict(patch, *cellIter);
-
-    // Compute action for cell
-    for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
-      const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
-      for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
-        const int iBlock = iBasis * spaceDim;
-        const double valI = wt*basis[iQ+iBasis];
-        for (int jBasis=0; jBasis < numBasis; ++jBasis) {
-          const int jBlock = jBasis * spaceDim;
-          const double val = valI * basis[iQ+jBasis];
-          for (int iDim=0; iDim < spaceDim; ++iDim)
-            _cellVector[iBlock+iDim] += val * fieldInCell[jBlock+iDim];
-        } // for
-      } // for
-    } // for
-    PetscErrorCode err = 
-      PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+2*spaceDim))));
-    if (err)
-      throw std::runtime_error("Logging PETSc flops failed.");
-    
-    // Assemble cell contribution into field
-    fieldOut->updateAdd(patch, *cellIter, _cellVector);
-  } // for
-} // integrateAction
-
-// ----------------------------------------------------------------------
-// Compute matrix associated with operator.
-void
-pylith::feassemble::IntegratorInertia::integrate(
-			     PetscMat* mat,
-			     const ALE::Obj<real_section_type>& fieldIn,
-			     const ALE::Obj<real_section_type>& coordinates)
-{ // integrate
-  assert(0 != mat);
-  assert(0 != _quadrature);
-  PetscErrorCode err;
-
-  // Get information about section
-  const topology_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = coordinates->getTopology();
-  const ALE::Obj<topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
-  const topology_type::label_sequence::iterator cellsEnd = cells->end();
-  const ALE::Obj<ALE::Mesh::order_type>& globalOrder = 
-    ALE::New::NumberingFactory<topology_type>::singleton(
-       topology->debug())->getGlobalOrder(topology, patch, 
-					  fieldIn->getName(), 
-					  fieldIn->getAtlas());
-
-  // Setup symmetric, sparse matrix
-  // :TODO: This needs to be moved outside Integrator object, because
-  // integrator object will be specific to cell type and material type
-  int localSize  = globalOrder->getLocalSize();
-  int globalSize = globalOrder->getGlobalSize();
-  err = MatCreate(topology->comm(), mat);
-  err = MatSetSizes(*mat, localSize, localSize, globalSize, globalSize);
-  err = MatSetFromOptions(*mat);
-  err = preallocateMatrix(topology, fieldIn->getAtlas(), globalOrder, *mat);
-
-  // Allocate matrix for cell values (if necessary)
-  _initCellMatrix();
-
-  for (topology_type::label_sequence::iterator cellIter=cells->begin();
-       cellIter != cellsEnd;
-       ++cellIter) {
-    // Compute geometry information for current cell
-    _quadrature->computeGeometry(coordinates, *cellIter);
-
-    // Reset element matrix to zero
-    _resetCellMatrix();
-
-    // Get cell geometry information
-    const int numQuadPts = _quadrature->numQuadPts();
-    const double* basis = _quadrature->basis();
-    const double* quadPts = _quadrature->quadPts();
-    const double* quadWts = _quadrature->quadWts();
-    const double* jacobianDet = _quadrature->jacobianDet();
-    const int numBasis = _quadrature->numCorners();
-    const int spaceDim = _quadrature->spaceDim();
-
-    // :TODO: Get mass density at quadrature points from material database
-    // For now, hardwire mass density
-    const double density = 1.0;
-
-    // Integrate cell
-    for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
-      const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density;
-      for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
-	const int iBlock = iBasis * spaceDim;
-	const double valI = wt*basis[iQ+iBasis];
-	for (int jBasis=0; jBasis < numBasis; ++jBasis) {
-	  const int jBlock = jBasis * spaceDim;
-	  const double val = valI * basis[iQ+jBasis];
-	  for (int iDim=0; iDim < spaceDim; ++iDim)
-	    _cellMatrix[(iBlock+iDim)*(numBasis*spaceDim)+jBlock+iDim] += val;
-	} // for
-      } // for
-    } // for
-    err = PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+spaceDim))));
-    if (err)
-      throw std::runtime_error("Logging PETSc flops failed.");
-    
-    // Assemble cell contribution into sparse matrix
-    err = updateOperator(*mat, fieldIn, globalOrder, *cellIter, _cellMatrix, 
-			 ADD_VALUES);
-  } // for
-} // integrate
-
-// ----------------------------------------------------------------------
-// Compute lumped matrix associated with operator.
-void
-pylith::feassemble::IntegratorInertia::integrateLumped(
-			     const ALE::Obj<real_section_type>& fieldOut,
-			     const ALE::Obj<real_section_type>& coordinates)
-{ // integrateLumped
-  assert(0 != _quadrature);
-
-  // Get information about section
-  const topology_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = coordinates->getTopology();
-  const ALE::Obj<topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
-  const topology_type::label_sequence::iterator cellsEnd = cells->end();
-
-  // Allocate matrix for cell values (if necessary)
-  _initCellVector();
-
-  for (topology_type::label_sequence::iterator cellIter=cells->begin();
-       cellIter != cellsEnd;
-       ++cellIter) {
-    // Compute geometry information for current cell
-    _quadrature->computeGeometry(coordinates, *cellIter);
-
-    // Reset element matrix to zero
-    _resetCellVector();
-
-    // Get cell geometry information
-    const int numQuadPts = _quadrature->numQuadPts();
-    const double* basis = _quadrature->basis();
-    const double* quadPts = _quadrature->quadPts();
-    const double* quadWts = _quadrature->quadWts();
-    const double* jacobianDet = _quadrature->jacobianDet();
-    const int numBasis = _quadrature->numCorners();
-    const int spaceDim = _quadrature->spaceDim();
-
-    // :TODO: Get mass density at quadrature points from material database
-    // For now, hardwire mass density
-    const double density = 1.0;
-
-    // Compute lumped mass matrix for cell
-    for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
-      const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density;
-      for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
-	const int iBlock = iBasis * spaceDim;
-	const double valI = wt*basis[iQ+iBasis];
-	for (int jBasis=0; jBasis < numBasis; ++jBasis) {
-	  const int jBlock = jBasis * spaceDim;
-	  const double val = valI*basis[iQ+jBasis];
-	  for (int iDim=0; iDim < spaceDim; ++iDim)
-	    _cellVector[iBlock+iDim] += val;
-	} // for
-      } // for
-    } // for
-
-    PetscErrorCode err = 
-      PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+spaceDim))));
-    if (err)
-      throw std::runtime_error("Logging PETSc flops failed.");
-    
-    // Assemble cell contribution into field
-    fieldOut->updateAdd(patch, *cellIter, _cellVector);
-  } // for
-} // integrateLumped
-
-// ----------------------------------------------------------------------
-// Initialize, get material property parameters from database.
-void
-pylith::feassemble::IntegratorInertia::initialize(
-				     ALE::Obj<ALE::Mesh>& mesh,
-				     spatialdata::geocoords::CoordSys* cs,
-				     spatialdata::spatialdb::SpatialDB* db)
-{ // initialize
-  assert(0 != cs);
-  assert(0 != db);
-
-  typedef ALE::Mesh::real_section_type real_section_type;
-  typedef ALE::Mesh::topology_type topology_type;
-
-  // Create density section
-  const int numQuadPts = _quadrature->numQuadPts();
-  const ALE::Mesh::int_section_type::patch_type patch = 0;
-  _density = mesh->getRealSection("density");
-  const int fiberDim = numQuadPts; // number of values in field per cell
-  _density->setName("density");
-  _density->setFiberDimensionByDepth(patch, 0, fiberDim);
-  _density->allocate();
-
-  // Open database
-  db->open();
-  const int numVals = 1;
-  const char* names[numVals];
-  names[0] = "density";
-  db->queryVals(names, numVals);
-  
-  const ALE::Obj<real_section_type>& coordinates = 
-    mesh->getRealSection("coordinates");
-  const ALE::Obj<topology_type>& topology = coordinates->getTopology();
-  const ALE::Obj<topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
-  const topology_type::label_sequence::iterator cellsEnd = cells->end();
-
-  // Loop over cells
-  double* cellDensity = (numQuadPts > 0) ? new double[numQuadPts] : 0;
-  for (topology_type::label_sequence::iterator cellIter=cells->begin();
-       cellIter != cellsEnd;
-       ++cellIter) {
-    // Compute geometry information for current cell
-    _quadrature->computeGeometry(coordinates, *cellIter);
-
-    const double* quadPts = _quadrature->quadPts();
-    const int spaceDim = _quadrature->spaceDim();
-
-    // Loop over quadrature points in cell and query database
-    for (int iQuadPt=0, index=0; 
-	 iQuadPt < numQuadPts; 
-	 ++iQuadPt, index+=spaceDim)
-      // account for differences in spaceDim
-      const int err = db->query(&cellDensity[iQuadPt], numVals, 
-				&quadPts[index], spaceDim, cs);
-    // Assemble cell contribution into field
-    _density->updateAdd(patch, *cellIter, cellDensity);
-  } // for
-  delete[] cellDensity; cellDensity = 0;
-
-  // Close database
-  db->close();
-} // initialize
-
-
-// End of file 

Deleted: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -1,127 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-//                           Brad T. Aagaard
-//                        U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-/**
- * @file pylith/feassemble/IntegratorInertia.hh
- *
- * @brief Integrate inertial term for 3-D finite elements.
- *
- * This object performs integrals over the domain of a finite-element
- * associated with translational inertia.
- *   - Integrate action over cell
- *     \f[
- *       \int_{V^e} \rho N^p \sum_q N^q u_i^q \, dV
- *     \f]
- *   - Integrate to form matrix
- *     \f[
- *       \int_{V^e} (\rho N^q N^q)_i \, dV
- *     \f]
- *   - Integrate and lump to form lumped matrix (field)
- *
- * See governing equations section of user manual for more
- * information.
- */
-
-#if !defined(pylith_feassemble_integratorinertia_hh)
-#define pylith_feassemble_integratorinertia_hh
-
-#include "Integrator.hh"
-
-namespace pylith {
-  namespace feassemble {
-    class IntegratorInertia;
-    class TestIntegratorInertia;
-  } // feassemble
-} // pylith
-
-class pylith::feassemble::IntegratorInertia : public Integrator
-{ // Integrator1D
-  friend class TestIntegratorInertia; // unit testing
-
-// PUBLIC MEMBERS ///////////////////////////////////////////////////////
-public :
-
-  /// Constructor
-  IntegratorInertia(void);
-
-  /// Destructor
-  ~IntegratorInertia(void);
-
-  /// Create a copy of this object.
-  Integrator* clone(void) const;
-
-  /** Integrate inertial term for 3-D finite elements.
-   *
-   * @param fieldOut Output field
-   * @param fieldIn Input field
-   * @param coordinates Field of cell vertex coordinates
-   */
-  void integrateAction(const ALE::Obj<ALE::Mesh::real_section_type>& fieldOut,
-		       const ALE::Obj<real_section_type>& fieldIn,
-		       const ALE::Obj<real_section_type>& coordinates);
-
-  /** Compute matrix associated with operator.
-   *
-   * @param mat Sparse matrix
-   * @param coordinates Field of cell vertex coordinates
-   */
-  void integrate(PetscMat* mat,
-                 const ALE::Obj<real_section_type>& fieldIn,
-                 const ALE::Obj<real_section_type>& coordinates);
-
-  /** Compute lumped mass matrix.
-   *
-   * @param fieldOut Lumped mass matrix as field
-   * @param coordinates Field of cell vertex coordinates
-   */
-  void integrateLumped(const ALE::Obj<real_section_type>& fieldOut,
-		       const ALE::Obj<real_section_type>& coordinates);
-
-  /** Initialize, get material property parameters from database.
-   *
-   * @param mesh PETSc mesh
-   * @param cs Pointer to coordinate system of vertices
-   * @param db Pointer to spatial database with material property parameters
-   */
-  void initialize(ALE::Obj<ALE::Mesh>& mesh,
-		  spatialdata::geocoords::CoordSys* cs,
-		  spatialdata::spatialdb::SpatialDB* db);
-
-// PROTECTED METHODS ////////////////////////////////////////////////////
-protected :
-
-  /** Copy constructor.
-   *
-   * @param i Integrator to copy
-   */
-  IntegratorInertia(const IntegratorInertia& i);
-
-// PRIVATE METHODS //////////////////////////////////////////////////////
-private :
-
-  /// Not implemented
-  const IntegratorInertia& operator=(const IntegratorInertia&);
-
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
-
-  /// Field with density at quadrature points
-  /// Fiber dimension is number of quadrature points per cell
-  ALE::Obj<real_section_type> _density;
-
-}; // IntegratorInertia
-
-#include "IntegratorInertia.icc" // inline methods
-
-#endif // pylith_feassemble_integratorinertia_hh
-
-// End of file 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/ParameterManager.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ParameterManager.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ParameterManager.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -46,13 +46,13 @@
   } // if
   
   ALE::Obj<real_section_type> parameter = 
-    new real_section_type(_mesh->getTopology());
+    new real_section_type(_mesh->comm(), _mesh->debug());
   _real[name] = parameter;
 } // addReal
 
 // ----------------------------------------------------------------------
 // Get parameter.
-const ALE::Obj<ALE::Mesh::real_section_type>&
+const ALE::Obj<pylith::feassemble::ParameterManager::Mesh::real_section_type>&
 pylith::feassemble::ParameterManager::getReal(const char* name)
 { // getReal
   map_real_type::const_iterator iter = _real.find(name);

Modified: short/3D/PyLith/trunk/libsrc/feassemble/ParameterManager.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ParameterManager.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ParameterManager.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -42,9 +42,8 @@
 // PUBLIC TYPEDEFS //////////////////////////////////////////////////////
 public :
 
-  typedef ALE::Mesh Mesh;
-  typedef Mesh::topology_type topology_type;
-  typedef topology_type::point_type point_type;
+  typedef ALE::Field::Mesh        Mesh;
+  typedef Mesh::point_type        point_type;
   typedef Mesh::real_section_type real_section_type;
 
 // PUBLIC MEMBERS ///////////////////////////////////////////////////////
@@ -87,7 +86,7 @@
 private :
 
   /// PETSc mesh associated with fields
-  const ALE::Obj<ALE::Mesh>& _mesh;
+  const ALE::Obj<Mesh>& _mesh;
 
   /// Map for parameters stored as real fields
   map_real_type _real;

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -46,8 +46,7 @@
 // PUBLIC TYPEDEFS //////////////////////////////////////////////////////
 public :
 
-  typedef ALE::Mesh Mesh;
-  typedef Mesh::topology_type topology_type;
+  typedef ALE::Field::Mesh        Mesh;
   typedef Mesh::real_section_type real_section_type;
 
 // PUBLIC METHODS ///////////////////////////////////////////////////////
@@ -186,8 +185,9 @@
    * @param cell Finite-element cell
    */
   virtual 
-  void computeGeometry(const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell) = 0;
+  void computeGeometry(const ALE::Obj<Mesh>& mesh,
+               const ALE::Obj<real_section_type>& coordinates,
+		       const Mesh::point_type& cell) = 0;
 
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -39,8 +39,9 @@
 // Compute geometric quantities for a cell.
 void
 pylith::feassemble::Quadrature1D::computeGeometry(
+		       const ALE::Obj<Mesh>& mesh,
 		       const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell)
+		       const Mesh::point_type& cell)
 { // computeGeometry
   assert(1 == _cellDim);
   assert(1 == _spaceDim);
@@ -54,9 +55,8 @@
   _resetGeometry();
 
   // Get coordinates of cell's vertices
-  const topology_type::patch_type patch  = 0;
   const real_section_type::value_type* vertCoords = 
-    coordinates->restrict(patch, cell);
+    mesh->restrict(coordinates, cell);
   //assert(1 == coordinates.GetFiberDimensionByDepth(patch, 
   //*vertices->begin(), 0));
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -44,8 +44,9 @@
   /// Create a copy of this object.
   Quadrature* clone(void) const;
 
-  void computeGeometry(const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell);
+  void computeGeometry(const ALE::Obj<Mesh>& mesh,
+               const ALE::Obj<real_section_type>& coordinates,
+		       const Mesh::point_type& cell);
 
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -39,8 +39,9 @@
 // Compute geometric quantities for a cell.
 void
 pylith::feassemble::Quadrature1Din2D::computeGeometry(
+		       const ALE::Obj<Mesh>& mesh,
 		       const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell)
+		       const Mesh::point_type& cell)
 { // computeGeometry
   assert(1 == _cellDim);
   assert(2 == _spaceDim);
@@ -54,9 +55,8 @@
   _resetGeometry();
 
   // Get coordinates of cell's vertices
-  const topology_type::patch_type patch  = 0;
   const real_section_type::value_type* vertCoords = 
-    coordinates->restrict(patch, cell);
+    mesh->restrict(coordinates, cell);
   //assert(2 == coordinates.GetFiberDimensionByDepth(patch, 
   //*vertices->begin(), 0));
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -49,8 +49,9 @@
    * @param coordinates Section containing vertex coordinates
    * @param cell Finite-element cell
    */
-  void computeGeometry(const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell);
+  void computeGeometry(const ALE::Obj<Mesh>& mesh,
+               const ALE::Obj<real_section_type>& coordinates,
+		       const Mesh::point_type& cell);
 
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -39,8 +39,9 @@
 // Compute geometric quantities for a cell.
 void
 pylith::feassemble::Quadrature1Din3D::computeGeometry(
+		       const ALE::Obj<Mesh>& mesh,
 		       const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell)
+		       const Mesh::point_type& cell)
 { // computeGeometry
   assert(1 == _cellDim);
   assert(3 == _spaceDim);
@@ -54,9 +55,8 @@
   _resetGeometry();
 
   // Get coordinates of cell's vertices
-  const ALE::Mesh::topology_type::patch_type patch  = 0;
   const ALE::Mesh::real_section_type::value_type* vertCoords = 
-    coordinates->restrict(patch, cell);
+    mesh->restrict(coordinates, cell);
   //assert(3 == coordinates.GetFiberDimensionByDepth(patch,
   //*vertices->begin(), 0));
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -49,8 +49,9 @@
    * @param coordinates Section containing vertex coordinates
    * @param cell Finite-element cell
    */
-  void computeGeometry(const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell);
+  void computeGeometry(const ALE::Obj<Mesh>& mesh,
+               const ALE::Obj<real_section_type>& coordinates,
+		       const Mesh::point_type& cell);
 
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -39,8 +39,9 @@
 // Compute geometric quantities for a cell.
 void
 pylith::feassemble::Quadrature2D::computeGeometry(
+		       const ALE::Obj<Mesh>& mesh,
 		       const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell)
+		       const Mesh::point_type& cell)
 { // computeGeometry
   assert(2 == _cellDim);
   assert(2 == _spaceDim);
@@ -54,9 +55,8 @@
   _resetGeometry();
 
   // Get coordinates of cell's vertices
-  const topology_type::patch_type patch  = 0;
   const real_section_type::value_type* vertCoords = 
-    coordinates->restrict(patch, cell);
+    mesh->restrict(coordinates, cell);
   //assert(2 == coordinates.GetFiberDimensionByDepth(patch,
   //*vertices->begin(), 0));
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -49,8 +49,9 @@
    * @param coordinates Section containing vertex coordinates
    * @param cell Finite-element cell
    */
-  void computeGeometry(const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell);
+  void computeGeometry(const ALE::Obj<Mesh>& mesh,
+               const ALE::Obj<real_section_type>& coordinates,
+		       const Mesh::point_type& cell);
 
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -42,8 +42,9 @@
 // Compute geometric quantities for a cell.
 void
 pylith::feassemble::Quadrature2Din3D::computeGeometry(
+		       const ALE::Obj<Mesh>& mesh,
 		       const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell)
+		       const Mesh::point_type& cell)
 { // computeGeometry
   assert(2 == _cellDim);
   assert(3 == _spaceDim);
@@ -57,9 +58,8 @@
   _resetGeometry();
 
   // Get coordinates of cell's vertices
-  const topology_type::patch_type patch  = 0;
   const real_section_type::value_type* vertCoords = 
-    coordinates->restrict(patch, cell);
+    mesh->restrict(coordinates, cell);
   //assert(3 == coordinates.GetFiberDimensionByDepth(patch,
   //*vertices->begin(), 0));
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -49,8 +49,9 @@
    * @param coordinates Section containing vertex coordinates
    * @param cell Finite-element cell
    */
-  void computeGeometry(const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell);
+  void computeGeometry(const ALE::Obj<Mesh>& mesh,
+               const ALE::Obj<real_section_type>& coordinates,
+		       const Mesh::point_type& cell);
 
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -39,8 +39,9 @@
 // Compute geometric quantities for a cell.
 void
 pylith::feassemble::Quadrature3D::computeGeometry(
+		       const ALE::Obj<Mesh>& mesh,
 		       const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell)
+		       const Mesh::point_type& cell)
 { // computeGeometry
   assert(3 == _cellDim);
   assert(3 == _spaceDim);
@@ -54,9 +55,8 @@
   _resetGeometry();
 
   // Get coordinates of cell's vertices
-  const topology_type::patch_type patch = 0;
   const real_section_type::value_type* vertCoords = 
-    coordinates->restrict(patch, cell);
+    mesh->restrict(coordinates, cell);
   //assert(3 == coordinates->getFiberDimensionByDepth(patch,
   //*vertices->begin(), 0));
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -49,8 +49,9 @@
    * @param coordinates Section containing vertex coordinates
    * @param cell Finite-element cell
    */
-  void computeGeometry(const ALE::Obj<real_section_type>& coordinates,
-		       const topology_type::point_type& cell);
+  void computeGeometry(const ALE::Obj<Mesh>& mesh,
+               const ALE::Obj<real_section_type>& coordinates,
+		       const Mesh::point_type& cell);
 
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -49,8 +49,7 @@
 // Compute physical properties of cell at quadrature points.
 void
 pylith::materials::ElasticMaterial::calcProperties(
-				     const topology_type::point_type& cell,
-				     const topology_type::patch_type& patch,
+				     const Mesh::point_type& cell,
 				     const int numQuadPts)
 { // calcProperties
   _initCellData(numQuadPts);
@@ -64,9 +63,9 @@
     const ALE::Obj<real_section_type> parameter = 
       _parameters->getReal(paramNames[iParam]);
 
-    assert(numQuadPts == parameter->getFiberDimension(patch, cell));
+    assert(numQuadPts == parameter->getFiberDimension(cell));
     const real_section_type::value_type* parameterCell =
-      parameter->restrict(patch, cell);
+      parameter->restrictPoint(cell);
     for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt)
       paramsCell[iQuadPt*numParams+iParam] = parameterCell[iQuadPt];
   } // for

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -41,8 +41,7 @@
   // PUBLIC TYPEDEFS ////////////////////////////////////////////////////
 public :
 
-  typedef ALE::Mesh Mesh;
-  typedef Mesh::topology_type topology_type;
+  typedef ALE::Field::Mesh        Mesh;
   typedef Mesh::real_section_type real_section_type;
 
   // PUBLIC METHODS /////////////////////////////////////////////////////
@@ -110,8 +109,7 @@
    * @param patch Finite-element patch
    * @param numQuadPts Number of quadrature points (consistency check)
    */
-  void calcProperties(const topology_type::point_type& cell,
-		      const topology_type::patch_type& patch,
+  void calcProperties(const Mesh::point_type& cell,
 		      const int numQuadPts);
 
   // PROTECTED METHODS //////////////////////////////////////////////////

Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -58,7 +58,7 @@
 // ----------------------------------------------------------------------
 // Get physical property parameters from database.
 void
-pylith::materials::Material::initialize(const ALE::Obj<ALE::Mesh>& mesh,
+pylith::materials::Material::initialize(const ALE::Obj<ALE::Field::Mesh>& mesh,
 					const spatialdata::geocoords::CoordSys* cs,
 					pylith::feassemble::Quadrature* quadrature)
 { // initialize
@@ -66,17 +66,14 @@
   assert(0 != cs);
   assert(0 != quadrature);
 
-  typedef ALE::Mesh::real_section_type real_section_type;
-  typedef ALE::Mesh::topology_type topology_type;
+  typedef ALE::Field::Mesh::real_section_type real_section_type;
 
   // Get cells associated with material
-  const ALE::Mesh::int_section_type::patch_type patch = 0;
   const ALE::Obj<real_section_type>& coordinates = 
     mesh->getRealSection("coordinates");
-  const ALE::Obj<topology_type>& topology = coordinates->getTopology();
-  const ALE::Obj<topology_type::label_sequence>& cells = 
-    topology->getLabelStratum(patch, "material-id", _id);
-  const topology_type::label_sequence::iterator cellsEnd = cells->end();
+  const ALE::Obj<ALE::Field::Mesh::label_sequence>& cells = 
+    mesh->getLabelStratum("material-id", _id);
+  const ALE::Field::Mesh::label_sequence::iterator cellsEnd = cells->end();
 
   // Check to make sure we have cells
   if (0 == cells->size()) {
@@ -99,8 +96,8 @@
   for (int iParam=0; iParam < numParams; ++iParam) {
     _parameters->addReal(paramNames[iParam]);
     paramSections[iParam] = _parameters->getReal(paramNames[iParam]);
-    paramSections[iParam]->setFiberDimension(patch, cells, fiberDim);
-    paramSections[iParam]->allocate();
+    paramSections[iParam]->setFiberDimension(cells, fiberDim);
+    mesh->allocate(paramSections[iParam]);
   } // for
 
   // Setup database for querying
@@ -114,11 +111,11 @@
   double** cellData = (numParams > 0) ? new double*[numParams] : 0;
   for (int iParam = 0; iParam < numParams; ++iParam)
     cellData[iParam] = (numQuadPts > 0) ? new double[numQuadPts] : 0;
-  for (topology_type::label_sequence::iterator cellIter=cells->begin();
+  for (ALE::Field::Mesh::label_sequence::iterator cellIter=cells->begin();
        cellIter != cellsEnd;
        ++cellIter) {
     // Compute geometry information for current cell
-    quadrature->computeGeometry(coordinates, *cellIter);
+    quadrature->computeGeometry(mesh, coordinates, *cellIter);
 
     const double* quadPts = quadrature->quadPts();
     const int spaceDim = quadrature->spaceDim();
@@ -154,7 +151,7 @@
     } // for
     // Assemble cell contribution into fields
     for (int iParam=0; iParam < numParams; ++iParam)
-      paramSections[iParam]->updateAdd(patch, *cellIter, cellData[iParam]);
+      mesh->updateAdd(paramSections[iParam], *cellIter, cellData[iParam]);
   } // for
   for (int iParam=0; iParam < numParams; ++iParam) {
     delete[] cellData[iParam]; cellData[iParam] = 0;

Modified: short/3D/PyLith/trunk/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -47,7 +47,9 @@
 
 /// Namespace for spatialdata package
 namespace ALE {
-  class Mesh;
+  namespace Field {
+    class Mesh;
+  }
   template<class T> class Obj;
 } // ALE
 
@@ -103,7 +105,7 @@
    * @param cs Coordinate system associated with mesh
    * @param quadrature Quadrature for finite-element integration
    */
-  void initialize(const ALE::Obj<ALE::Mesh>& mesh,
+  void initialize(const ALE::Obj<ALE::Field::Mesh>& mesh,
 		  const spatialdata::geocoords::CoordSys* cs,
 		  pylith::feassemble::Quadrature* quadrature);
   

Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -84,20 +84,15 @@
   ALE::Obj<Mesh>& mesh = *_mesh;
   
   ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
-  ALE::Obj<topology_type> topology = new topology_type(mesh->comm());
 
-  ALE::New::SieveBuilder<sieve_type>::buildTopology(sieve, meshDim, 
+  ALE::New::SieveBuilder<Mesh>::buildTopology(sieve, meshDim, 
 						    numCells, 
 						    const_cast<int*>(cells), 
 						    numVertices, 
 						    _interpolate, numCorners);
-  sieve->stratify();
-  topology->setPatch(0, sieve);
-  topology->stratify();
-  mesh->setTopology(topology);
-  ALE::New::SieveBuilder<sieve_type>::buildCoordinates(
-		      mesh->getRealSection("coordinates"), 
-		      spaceDim, coordinates);
+  mesh->setSieve(sieve);
+  mesh->stratify();
+  ALE::New::SieveBuilder<Mesh>::buildCoordinatesNew(mesh, spaceDim, coordinates);
 } // _buildMesh
 
 // ----------------------------------------------------------------------
@@ -110,17 +105,13 @@
   assert(0 != _mesh);
   ALE::Obj<Mesh>& mesh = *_mesh;
 
-  const Mesh::real_section_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = mesh->getTopology();
-
-  const ALE::Obj<Mesh::topology_type::label_sequence>& vertices = 
-    topology->depthStratum(patch, 0);
+  const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
   const ALE::Obj<Mesh::real_section_type>& coordsField =
     mesh->getRealSection("coordinates");
 
   const int numVertices = vertices->size();
   const int spaceDim = 
-    coordsField->getFiberDimension(patch, *vertices->begin());
+    coordsField->getFiberDimension(*vertices->begin());
 
   double* coordinates = 0;
   const int size = numVertices * spaceDim;
@@ -128,12 +119,11 @@
     coordinates = new double[size];
 
     int i = 0;
-    for(Mesh::topology_type::label_sequence::iterator v_iter = 
-	  vertices->begin();
-	v_iter != vertices->end();
-	++v_iter) {
+    for(Mesh::label_sequence::iterator v_iter = vertices->begin();
+        v_iter != vertices->end();
+	    ++v_iter) {
       const Mesh::real_section_type::value_type *vertexCoords = 
-	coordsField->restrict(patch, *v_iter);
+	coordsField->restrictPoint(*v_iter);
       for (int iDim=0; iDim < spaceDim; ++iDim)
 	coordinates[i++] = vertexCoords[iDim];
     } // for
@@ -158,20 +148,16 @@
   assert(0 != _mesh);
   ALE::Obj<Mesh>& mesh = *_mesh;
 
-  const topology_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = mesh->getTopology();
+  const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
+  const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
 
-  const ALE::Obj<sieve_type>& sieve = topology->getPatch(patch);
-  const ALE::Obj<Mesh::topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
-
   const int meshDim = mesh->getDimension();
   const int numCells = cells->size();
   const int numCorners = sieve->nCone(*cells->begin(), 
-				      topology->depth())->size();
+				      mesh->depth())->size();
 
   const ALE::Obj<Mesh::numbering_type>& vNumbering = 
-    mesh->getFactory()->getLocalNumbering(topology, patch, 0);
+    mesh->getFactory()->getLocalNumbering(mesh, 0);
 
   int* cellsArray = 0;
   const int size = numCells * numCorners;
@@ -180,7 +166,7 @@
     
     const int offset = (useIndexZero()) ? 0 : 1;
     int i = 0;
-    for(Mesh::topology_type::label_sequence::iterator e_iter = cells->begin();
+    for(Mesh::label_sequence::iterator e_iter = cells->begin();
 	e_iter != cells->end();
 	++e_iter) {
       const ALE::Obj<sieve_type::traits::coneSequence>& cone = 
@@ -211,10 +197,7 @@
   assert(0 != _mesh);
   ALE::Obj<Mesh>& mesh = *_mesh;
   
-  const topology_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = mesh->getTopology();
-  const ALE::Obj<Mesh::topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
+  const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
 
   if (cells->size() != numCells) {
     std::ostringstream msg;
@@ -224,14 +207,13 @@
     throw std::runtime_error(msg.str());
   } // if
 
-  const ALE::Obj<patch_label_type>& labelMaterials = 
-    topology->createLabel(patch, "material-id");
+  const ALE::Obj<label_type>& labelMaterials = mesh->createLabel("material-id");
   
   int i = 0;
-  for(Mesh::topology_type::label_sequence::iterator e_iter = cells->begin();
+  for(Mesh::label_sequence::iterator e_iter = cells->begin();
       e_iter != cells->end();
       ++e_iter)
-    topology->setValue(labelMaterials, *e_iter, materialIds[i++]);
+    mesh->setValue(labelMaterials, *e_iter, materialIds[i++]);
 } // _setMaterials
 
 // ----------------------------------------------------------------------
@@ -243,10 +225,7 @@
   assert(0 != _mesh);
   ALE::Obj<Mesh>& mesh = *_mesh;
 
-  const topology_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = mesh->getTopology();
-  const ALE::Obj<Mesh::topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
+  const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
   const int numCells = cells->size();
 
   int* materialsArray = 0;
@@ -254,16 +233,15 @@
   if (0 != pMaterialIds && size > 0) {
     materialsArray = new int[size];
   
-    const ALE::Obj<patch_label_type>& labelMaterials = 
-      topology->getLabel(patch, "material-id");
+    const ALE::Obj<label_type>& labelMaterials = mesh->getLabel("material-id");
     const int idDefault = 0;
 
     int i = 0;
-    for(Mesh::topology_type::label_sequence::iterator e_iter = cells->begin();
+    for(Mesh::label_sequence::iterator e_iter = cells->begin();
 	e_iter != cells->end();
 	++e_iter)
       materialsArray[i++] = 
-	topology->getValue(labelMaterials, *e_iter, idDefault);
+	mesh->getValue(labelMaterials, *e_iter, idDefault);
   } // if  
 
   if (0 != pMaterialIds)

Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -25,10 +25,9 @@
 { // MeshIO
 // PUBLIC TYPEDEFS //////////////////////////////////////////////////////
 public :
-  typedef ALE::Mesh Mesh;
-  typedef ALE::Mesh::sieve_type sieve_type;
-  typedef ALE::Mesh::topology_type topology_type;
-  typedef ALE::Sifter<int, sieve_type::point_type, int> patch_label_type;
+  typedef ALE::Field::Mesh Mesh;
+  typedef Mesh::sieve_type sieve_type;
+  typedef Mesh::label_type label_type;
   
 // PUBLIC MEMBERS ///////////////////////////////////////////////////////
 public :

Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src	2007-03-30 23:49:16 UTC (rev 6476)
@@ -441,21 +441,25 @@
     return
 
 
-  def integrateConstant(self, fieldOut, fieldInT, fieldInTmdt, coords):
+  def integrateConstant(self, mesh, fieldOut, fieldInT, fieldInTmdt, coords):
     """
     Integrate constant term (b) for dynamic elasticity term for 3-D
     finite elements.
     """
     # create shim for method 'integrateConstant'
-    #embed{ void IntegratorExplicit_integrateConstant(void* objVptr, void* fieldOutVptr, void* fieldInTVptr, void* fieldInTmdtVptr, void* coordsVptr)
-    typedef ALE::Mesh::real_section_type real_section_type;
+    #embed{ void IntegratorExplicit_integrateConstant(void* objVptr, void* meshVptr, void* fieldOutVptr, void* fieldInTVptr, void* fieldInTmdtVptr, void* coordsVptr)
+    typedef ALE::Field::Mesh Mesh;
+    typedef ALE::Field::Mesh::real_section_type real_section_type;
 
     try {
       assert(0 != objVptr);
+      assert(0 != meshVptr);
       assert(0 != fieldOutVptr);
       assert(0 != fieldInTVptr);
       assert(0 != fieldInTmdtVptr);
       assert(0 != coordsVptr);
+      ALE::Obj<Mesh>* mesh =
+        (ALE::Obj<Mesh>*) meshVptr;
       ALE::Obj<real_section_type>* fieldOut =
         (ALE::Obj<real_section_type>*) fieldOutVptr;
       ALE::Obj<real_section_type>* fieldInT =
@@ -465,7 +469,7 @@
       ALE::Obj<real_section_type>* coords =
         (ALE::Obj<real_section_type>*) coordsVptr;
       ((pylith::feassemble::IntegratorExplicit*) objVptr)->integrateConstant(
-                *fieldOut, *fieldInT, *fieldInTmdt, *coords);
+                *mesh, *fieldOut, *fieldInT, *fieldInTmdt, *coords);
     } catch (const std::exception& err) {
       PyErr_SetString(PyExc_RuntimeError,
                       const_cast<char*>(err.what()));
@@ -477,21 +481,23 @@
                       "Caught unknown C++ exception.");
     } // try/catch
     #}embed
+    cdef void* meshVptr
     cdef void* fieldOutVptr
     cdef void* fieldInTVptr
     cdef void* fieldInTmdtVptr
     cdef void* coordsVptr
+    meshVptr = PyCObject_AsVoidPtr(mesh)
     fieldOutVptr = PyCObject_AsVoidPtr(fieldOut)
     fieldInTVptr = PyCObject_AsVoidPtr(fieldInT)
     fieldInTmdtVptr = PyCObject_AsVoidPtr(fieldInTmdt)
     coordsVptr = PyCObject_AsVoidPtr(coords)
-    IntegratorExplicit_integrateConstant(self.thisptr,
+    IntegratorExplicit_integrateConstant(self.thisptr, meshVptr,
                                          fieldOutVptr, fieldInTVptr,
                                          fieldInTmdtVptr, coordsVptr)
     return
 
 
-  def integrateJacobian(self, mat, fieldIn, coords):
+  def integrateJacobian(self, mat, mesh, fieldIn, coords):
     """
     Compute matrix (A) associated with operator.
     """

Modified: short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src	2007-03-30 23:49:16 UTC (rev 6476)
@@ -71,7 +71,7 @@
     try {
       assert(0 != pObj);
       assert(0 != pMeshObj);
-      ALE::Obj<ALE::Mesh>* pMesh = (ALE::Obj<ALE::Mesh>*) pMeshObj;
+      ALE::Obj<ALE::Field::Mesh>* pMesh = (ALE::Obj<ALE::Field::Mesh>*) pMeshObj;
       ((pylith::meshio::MeshIO*) pObj)->read(pMesh);
     } catch (const std::exception& err) {
       PyErr_SetString(PyExc_RuntimeError,
@@ -95,7 +95,7 @@
     # create shim for method 'write'
     #embed{ void MeshIO_write(void* pObj, void* pMeshObj)
     try {
-      ALE::Obj<ALE::Mesh>* pMesh = (ALE::Obj<ALE::Mesh>*) pMeshObj;
+      ALE::Obj<ALE::Field::Mesh>* pMesh = (ALE::Obj<ALE::Field::Mesh>*) pMeshObj;
       ((pylith::meshio::MeshIO*) pObj)->write(pMesh);
     } catch (const std::exception& err) {
       PyErr_SetString(PyExc_RuntimeError,

Modified: short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src	2007-03-30 23:49:16 UTC (rev 6476)
@@ -37,7 +37,7 @@
   Destroy Petsc Mesh.
   """
   #embed{ void MeshPtr_destructor_cpp(void* pObj)
-  ALE::Obj<ALE::Mesh>* pMesh = (ALE::Obj<ALE::Mesh>*) pObj;
+  ALE::Obj<ALE::Field::Mesh>* pMesh = (ALE::Obj<ALE::Field::Mesh>*) pObj;
   delete pMesh;
   #}embed
   MeshPtr_destructor_cpp(obj)
@@ -58,7 +58,7 @@
     #embed{ void* MeshPtr_constructor()
     void* result = 0;
     try {
-      result = (void*)(new ALE::Obj<ALE::Mesh>);
+      result = (void*)(new ALE::Obj<ALE::Field::Mesh>);
     } catch (const std::exception& err) {
       PyErr_SetString(PyExc_RuntimeError,
                       const_cast<char*>(err.what()));
@@ -81,11 +81,11 @@
     """
     # create shim for getRealSection
     #embed{ void* Mesh_getRealSection(void* pObj, char* label)
-    typedef ALE::Mesh::real_section_type real_section_type;
+    typedef ALE::Field::Mesh::real_section_type real_section_type;
     
     void* result = 0;
     try {
-      ALE::Obj<ALE::Mesh>* mesh = (ALE::Obj<ALE::Mesh>*) pObj;
+      ALE::Obj<ALE::Field::Mesh>* mesh = (ALE::Obj<ALE::Field::Mesh>*) pObj;
       assert(0 != mesh);
       const ALE::Obj<real_section_type>& section =
         (*mesh)->getRealSection(label);
@@ -111,18 +111,17 @@
     """
     # create shim for createRealSection
     #embed{ void* Mesh_createRealSection(void* pObj, char* label, int fiberDim)
-    typedef ALE::Mesh::real_section_type real_section_type;
+    typedef ALE::Field::Mesh::real_section_type real_section_type;
     
     void* result = 0;
     try {
-      ALE::Obj<ALE::Mesh>* mesh = (ALE::Obj<ALE::Mesh>*) pObj;
+      ALE::Obj<ALE::Field::Mesh>* mesh = (ALE::Obj<ALE::Field::Mesh>*) pObj;
       assert(0 != mesh);
       const ALE::Obj<real_section_type>& section =
         (*mesh)->getRealSection(label);
       assert(!section.isNull());
-      const ALE::Mesh::topology_type::patch_type patch = 0;      
-      section->setFiberDimensionByDepth(patch, 0, fiberDim);
-      section->allocate();
+      section->setFiberDimension((*mesh)->depthStratum(0), fiberDim);
+      (*mesh)->allocate(section);
       result = (void*) &section;
     } catch (const std::exception& err) {
       PyErr_SetString(PyExc_RuntimeError,
@@ -151,14 +150,13 @@
   """
   # create shim for zero section
   #embed{ void* Section_zero(void* pObj)
-  typedef ALE::Mesh::real_section_type real_section_type;
+  typedef ALE::Field::Mesh::real_section_type real_section_type;
   
   try {
     ALE::Obj<real_section_type>* section =
       (ALE::Obj<real_section_type>*) pObj;
     assert(!section->isNull());
-    const ALE::Mesh::topology_type::patch_type patch = 0;      
-    (*section)->zero(patch);
+    (*section)->zero();
   } catch (const std::exception& err) {
     PyErr_SetString(PyExc_RuntimeError,
                     const_cast<char*>(err.what()));

Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -232,30 +232,24 @@
 		    cellDim, numCorners, numQuadPts, spaceDim);
 
   // Create mesh with test cell
-  typedef ALE::Mesh::topology_type topology_type;
-  typedef topology_type::sieve_type sieve_type;
-  ALE::Obj<ALE::Mesh> mesh = new ALE::Mesh(PETSC_COMM_WORLD, cellDim);
+  typedef ALE::Field::Mesh Mesh;
+  typedef ALE::Field::Mesh::sieve_type sieve_type;
+  ALE::Obj<Mesh> mesh = new Mesh(PETSC_COMM_WORLD, cellDim);
   ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
-  ALE::Obj<topology_type> topology = new topology_type(mesh->comm());
 
   const bool interpolate = false;
-  ALE::New::SieveBuilder<sieve_type>::buildTopology(sieve, cellDim, numCells,
+  ALE::New::SieveBuilder<Mesh>::buildTopology(sieve, cellDim, numCells,
 		     (int*) cells, numVertices, interpolate, numCorners);
-  sieve->stratify();
-  topology->setPatch(0, sieve);
-  topology->stratify();
-  mesh->setTopology(topology);
-  ALE::New::SieveBuilder<sieve_type>::buildCoordinates(
-		    mesh->getRealSection("coordinates"), spaceDim, vertCoords);
+  mesh->setSieve(sieve);
+  mesh->stratify();
+  ALE::New::SieveBuilder<Mesh>::buildCoordinatesNew(mesh, spaceDim, vertCoords);
   
   // Check values from computeGeometry()
-  const ALE::Mesh::topology_type::patch_type patch = 0;
-  const ALE::Obj<topology_type::label_sequence>& elements = 
-    topology->heightStratum(patch, 0);
-  const topology_type::label_sequence::iterator e_iter = elements->begin(); 
-  const ALE::Obj<ALE::Mesh::real_section_type>& coordinates = 
+  const ALE::Obj<Mesh::label_sequence>& elements = mesh->heightStratum(0);
+  const Mesh::label_sequence::iterator e_iter = elements->begin(); 
+  const ALE::Obj<Mesh::real_section_type>& coordinates = 
     mesh->getRealSection("coordinates");
-  pQuad->computeGeometry(coordinates, *e_iter);
+  pQuad->computeGeometry(mesh, coordinates, *e_iter);
 
   CPPUNIT_ASSERT(1 == numCells);
 

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -85,11 +85,11 @@
 void
 pylith::materials::TestElasticMaterial::testCalcProperties(void)
 { // testCalcProperties
-  typedef ALE::Mesh::topology_type topology_type;
-  typedef topology_type::sieve_type sieve_type;
-  typedef ALE::Mesh::real_section_type real_section_type;
+  typedef ALE::Field::Mesh Mesh;
+  typedef Mesh::sieve_type sieve_type;
+  typedef Mesh::real_section_type real_section_type;
 
-  ALE::Obj<ALE::Mesh> mesh;
+  ALE::Obj<Mesh> mesh;
   { // create mesh
     const int cellDim = 1;
     const int numCorners = 2;
@@ -101,28 +101,21 @@
     CPPUNIT_ASSERT(0 != vertCoords);
     CPPUNIT_ASSERT(0 != cells);
 
-    mesh = new ALE::Mesh(PETSC_COMM_WORLD, cellDim);
+    mesh = new Mesh(PETSC_COMM_WORLD, cellDim);
     ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
-    ALE::Obj<topology_type> topology = new topology_type(mesh->comm());
 
     const bool interpolate = false;
-    ALE::New::SieveBuilder<sieve_type>::buildTopology(sieve, cellDim, numCells,
+    ALE::New::SieveBuilder<Mesh>::buildTopology(sieve, cellDim, numCells,
 	       const_cast<int*>(cells), numVertices, interpolate, numCorners);
-    sieve->stratify();
-    topology->setPatch(0, sieve);
-    topology->stratify();
-    mesh->setTopology(topology);
-    ALE::New::SieveBuilder<sieve_type>::buildCoordinates(
-		   mesh->getRealSection("coordinates"), spaceDim, vertCoords);
+    mesh->setSieve(sieve);
+    mesh->stratify();
+    ALE::New::SieveBuilder<Mesh>::buildCoordinatesNew(mesh, spaceDim, vertCoords);
   } // create mesh
 
   // Get cells associated with material
-  const ALE::Mesh::int_section_type::patch_type patch = 0;
   const ALE::Obj<real_section_type>& coordinates = 
     mesh->getRealSection("coordinates");
-  const ALE::Obj<topology_type>& topology = coordinates->getTopology();
-  const ALE::Obj<topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
+  const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
 
   ElasticIsotropic3D material;
   ElasticIsotropic3DData data;
@@ -131,37 +124,37 @@
   const int numQuadPts = 2;
   const int fiberDim = numQuadPts; // number of values in field per cell
 
-  topology_type::label_sequence::iterator cellIter=cells->begin();
+  Mesh::label_sequence::iterator cellIter=cells->begin();
 
   material._parameters->addReal("density");
   const ALE::Obj<real_section_type>& parameterDensity = 
     material._parameters->getReal("density");
-  parameterDensity->setFiberDimension(patch, cells, fiberDim);
-  parameterDensity->allocate();
+  parameterDensity->setFiberDimension(cells, fiberDim);
+  mesh->allocate(parameterDensity);
   double cellData[numQuadPts];
   cellData[0] = data.parameterData[0];
   cellData[1] = data.parameterData[3];
-  parameterDensity->updateAdd(patch, *cellIter, cellData);
+  parameterDensity->updateAddPoint(*cellIter, cellData);
 
   material._parameters->addReal("mu");
   const ALE::Obj<real_section_type>& parameterMu = 
     material._parameters->getReal("mu");
-  parameterMu->setFiberDimension(patch, cells, fiberDim);
-  parameterMu->allocate();
+  parameterMu->setFiberDimension(cells, fiberDim);
+  mesh->allocate(parameterMu);
   cellData[0] = data.parameterData[1];
   cellData[1] = data.parameterData[4];
-  parameterMu->updateAdd(patch, *cellIter, cellData);
+  parameterMu->updateAddPoint(*cellIter, cellData);
 
   material._parameters->addReal("lambda");
   const ALE::Obj<real_section_type>& parameterLambda = 
     material._parameters->getReal("lambda");
-  parameterLambda->setFiberDimension(patch, cells, fiberDim);
-  parameterLambda->allocate();
+  parameterLambda->setFiberDimension(cells, fiberDim);
+  mesh->allocate(parameterLambda);
   cellData[0] = data.parameterData[2];
   cellData[1] = data.parameterData[5];
-  parameterLambda->updateAdd(patch, *cellIter, cellData);
+  parameterLambda->updateAddPoint(*cellIter, cellData);
 
-  material.calcProperties(*cellIter, patch, numQuadPts);
+  material.calcProperties(*cellIter, numQuadPts);
 
   const double tolerance = 1.0e-06;
 

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -75,12 +75,12 @@
 void
 pylith::materials::TestMaterial::testInitialize(void)
 { // testInitialize
-  typedef ALE::Mesh::topology_type topology_type;
-  typedef topology_type::sieve_type sieve_type;
-  typedef ALE::Mesh::real_section_type real_section_type;
-  typedef ALE::Sifter<int, sieve_type::point_type, int> patch_label_type;
+  typedef ALE::Field::Mesh Mesh;
+  typedef Mesh::sieve_type sieve_type;
+  typedef Mesh::label_type label_type;
+  typedef Mesh::real_section_type real_section_type;
 
-  ALE::Obj<ALE::Mesh> mesh;
+  ALE::Obj<Mesh> mesh;
   const int materialID = 24;
   { // create mesh
     const int cellDim = 1;
@@ -93,34 +93,26 @@
     CPPUNIT_ASSERT(0 != vertCoords);
     CPPUNIT_ASSERT(0 != cells);
 
-    mesh = new ALE::Mesh(PETSC_COMM_WORLD, cellDim);
+    mesh = new Mesh(PETSC_COMM_WORLD, cellDim);
     ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
-    ALE::Obj<topology_type> topology = new topology_type(mesh->comm());
 
     const bool interpolate = false;
-    ALE::New::SieveBuilder<sieve_type>::buildTopology(sieve, cellDim, numCells,
+    ALE::New::SieveBuilder<Mesh>::buildTopology(sieve, cellDim, numCells,
 	       const_cast<int*>(cells), numVertices, interpolate, numCorners);
-    sieve->stratify();
-    topology->setPatch(0, sieve);
-    topology->stratify();
-    mesh->setTopology(topology);
-    ALE::New::SieveBuilder<sieve_type>::buildCoordinates(
-		   mesh->getRealSection("coordinates"), spaceDim, vertCoords);
+    mesh->setSieve(sieve);
+    mesh->stratify();
+    ALE::New::SieveBuilder<Mesh>::buildCoordinatesNew(mesh, spaceDim, vertCoords);
 
   } // create mesh
 
   { // set material ids
-    const topology_type::patch_type patch = 0;
-    const ALE::Obj<topology_type>& topology = mesh->getTopology();
-    const ALE::Obj<ALE::Mesh::topology_type::label_sequence>& cells = 
-      topology->heightStratum(patch, 0);
-    const ALE::Obj<patch_label_type>& labelMaterials = 
-      topology->createLabel(patch, "material-id");
+    const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+    const ALE::Obj<label_type>& labelMaterials = mesh->createLabel("material-id");
     int i = 0;
-    for(ALE::Mesh::topology_type::label_sequence::iterator e_iter = cells->begin();
+    for(Mesh::label_sequence::iterator e_iter = cells->begin();
 	e_iter != cells->end();
 	++e_iter)
-      topology->setValue(labelMaterials, *e_iter, materialID);
+      mesh->setValue(labelMaterials, *e_iter, materialID);
   } // set material ids
 
   spatialdata::geocoords::CSCart cs;
@@ -166,18 +158,15 @@
   const double lambdaE[] = { lambdaA, lambdaB };
 
   // Get cells associated with material
-  const ALE::Mesh::int_section_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = mesh->getTopology();
-  const ALE::Obj<topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
+  const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
 
-  topology_type::label_sequence::iterator cellIter=cells->begin();
+  Mesh::label_sequence::iterator cellIter=cells->begin();
   const double tolerance = 1.0e-06;
 
   const ALE::Obj<real_section_type>& parameterDensity = 
     material._parameters->getReal("density");
   const real_section_type::value_type* densityCell = 
-    parameterDensity->restrict(patch, *cellIter);
+    parameterDensity->restrictPoint(*cellIter);
   CPPUNIT_ASSERT(0 != densityCell);
   for (int i=0; i < numQuadPts; ++i)
     CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, densityCell[i]/densityE[i], tolerance);
@@ -185,7 +174,7 @@
   const ALE::Obj<real_section_type>& parameterMu = 
     material._parameters->getReal("mu");
   const real_section_type::value_type* muCell = 
-    parameterMu->restrict(patch, *cellIter);
+    parameterMu->restrictPoint(*cellIter);
   CPPUNIT_ASSERT(0 != muCell);
   for (int i=0; i < numQuadPts; ++i)
     CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, muCell[i]/muE[i], tolerance);
@@ -193,7 +182,7 @@
   const ALE::Obj<real_section_type>& parameterLambda = 
     material._parameters->getReal("lambda");
   const real_section_type::value_type* lambdaCell = 
-    parameterLambda->restrict(patch, *cellIter);
+    parameterLambda->restrictPoint(*cellIter);
   CPPUNIT_ASSERT(0 != lambdaCell);
   for (int i=0; i < numQuadPts; ++i)
     CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, lambdaCell[i]/lambdaE[i], tolerance);

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc	2007-03-30 23:49:16 UTC (rev 6476)
@@ -20,12 +20,12 @@
 
 // ----------------------------------------------------------------------
 // Get simple mesh for testing I/O.
-ALE::Obj<ALE::Mesh>*
+ALE::Obj<ALE::Field::Mesh>*
 pylith::meshio::TestMeshIO::createMesh(const MeshData& data)
 { // createMesh
-  typedef ALE::Mesh::topology_type topology_type;
-  typedef topology_type::sieve_type sieve_type;
-  typedef ALE::Sifter<int, sieve_type::point_type, int> patch_label_type;
+  typedef ALE::Field::Mesh Mesh;
+  typedef Mesh::sieve_type sieve_type;
+  typedef Mesh::label_type label_type;
 
   // buildTopology() requires zero based index
   assert(true == data.useIndexZero);
@@ -42,35 +42,28 @@
   CPPUNIT_ASSERT(0 != cells);
   CPPUNIT_ASSERT(0 != materialIds);
 
-  ALE::Obj<ALE::Mesh>* meshHandle = new ALE::Obj<ALE::Mesh>;
-  *meshHandle = new ALE::Mesh(PETSC_COMM_WORLD, cellDim);
-  ALE::Obj<ALE::Mesh> mesh = *meshHandle;
+  ALE::Obj<Mesh>* meshHandle = new ALE::Obj<Mesh>;
+  *meshHandle = new Mesh(PETSC_COMM_WORLD, cellDim);
+  ALE::Obj<Mesh> mesh = *meshHandle;
   ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
-  ALE::Obj<topology_type> topology = new topology_type(mesh->comm());
 
   const bool interpolate = false;
-  ALE::New::SieveBuilder<sieve_type>::buildTopology(sieve, cellDim, numCells,
+  ALE::New::SieveBuilder<Mesh>::buildTopology(sieve, cellDim, numCells,
 	       const_cast<int*>(cells), numVertices, interpolate, numCorners);
-  sieve->stratify();
-  topology->setPatch(0, sieve);
-  topology->stratify();
-  mesh->setTopology(topology);
-  ALE::New::SieveBuilder<sieve_type>::buildCoordinates(
-		    mesh->getRealSection("coordinates"), spaceDim, vertCoords);
+  mesh->setSieve(sieve);
+  mesh->stratify();
+  ALE::New::SieveBuilder<Mesh>::buildCoordinatesNew(mesh, spaceDim, vertCoords);
 
-  const Mesh::real_section_type::patch_type patch = 0;
-  const ALE::Obj<Mesh::topology_type::label_sequence>& cellsMesh = 
-    topology->heightStratum(patch, 0);
+  const ALE::Obj<Mesh::label_sequence>& cellsMesh = mesh->heightStratum(0);
 
-  const ALE::Obj<patch_label_type>& labelMaterials = 
-    topology->createLabel(patch, "material-id");
+  const ALE::Obj<label_type>& labelMaterials = mesh->createLabel("material-id");
   
   int i = 0;
-  for(Mesh::topology_type::label_sequence::iterator e_iter = 
+  for(Mesh::label_sequence::iterator e_iter = 
 	cellsMesh->begin();
       e_iter != cellsMesh->end();
       ++e_iter)
-    topology->setValue(labelMaterials, *e_iter, materialIds[i++]);
+    mesh->setValue(labelMaterials, *e_iter, materialIds[i++]);
 
   return meshHandle;
 } // createMesh
@@ -78,58 +71,57 @@
 // ----------------------------------------------------------------------
 // Check values in mesh against data.
 void
-pylith::meshio::TestMeshIO::checkVals(const ALE::Obj<ALE::Mesh>& mesh,
+pylith::meshio::TestMeshIO::checkVals(const ALE::Obj<Mesh>& mesh,
 				      const MeshData& data)
 { // checkVals
-  typedef ALE::Sifter<int, sieve_type::point_type, int> patch_label_type;
+  typedef ALE::Field::Mesh::label_type label_type;
 
-  const Mesh::real_section_type::patch_type patch = 0;
-  const ALE::Obj<topology_type>& topology = mesh->getTopology();
-
   // Check mesh dimension
   CPPUNIT_ASSERT_EQUAL(data.cellDim, mesh->getDimension());
 
   // Check vertices
-  const ALE::Obj<Mesh::topology_type::label_sequence>& vertices = 
-    topology->depthStratum(patch, 0);
+  const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
   const ALE::Obj<Mesh::real_section_type>& coordsField =
     mesh->getRealSection("coordinates");
   const int numVertices = vertices->size();
   CPPUNIT_ASSERT_EQUAL(data.numVertices, numVertices);
   CPPUNIT_ASSERT_EQUAL(data.spaceDim, 
-		       coordsField->getFiberDimension(patch, 
-						      *vertices->begin()));
+		       coordsField->getFiberDimension(*vertices->begin()));
   int i = 0;
   const int spaceDim = data.spaceDim;
-  for(Mesh::topology_type::label_sequence::iterator v_iter = 
+  for(Mesh::label_sequence::iterator v_iter = 
 	vertices->begin();
       v_iter != vertices->end();
       ++v_iter) {
     const Mesh::real_section_type::value_type *vertexCoords = 
-      coordsField->restrict(patch, *v_iter);
+      coordsField->restrictPoint(*v_iter);
     const double tolerance = 1.0e-06;
     for (int iDim=0; iDim < spaceDim; ++iDim)
-      CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vertexCoords[iDim]/data.vertices[i++],
+      if (data.vertices[i] < 1.0) {
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(data.vertices[i++], vertexCoords[iDim],
 				   tolerance);
+      } else {
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vertexCoords[iDim]/data.vertices[i++],
+				   tolerance);
+      }
   } // for
 
   // check cells
-  const ALE::Obj<sieve_type>& sieve = topology->getPatch(patch);
-  const ALE::Obj<Mesh::topology_type::label_sequence>& cells = 
-    topology->heightStratum(patch, 0);
+  const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
+  const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
 
   const int numCells = cells->size();
   CPPUNIT_ASSERT_EQUAL(data.numCells, numCells);
   const int numCorners = sieve->nCone(*cells->begin(), 
-				      topology->depth())->size();
+				      mesh->depth())->size();
   CPPUNIT_ASSERT_EQUAL(data.numCorners, numCorners);
 
   const ALE::Obj<Mesh::numbering_type>& vNumbering = 
-    mesh->getFactory()->getLocalNumbering(topology, patch, 0);
+    mesh->getFactory()->getLocalNumbering(mesh, 0);
 
   const int offset = (data.useIndexZero) ? 0 : 1;
   i = 0;
-  for(Mesh::topology_type::label_sequence::iterator e_iter = cells->begin();
+  for(Mesh::label_sequence::iterator e_iter = cells->begin();
       e_iter != cells->end();
       ++e_iter) {
     const ALE::Obj<sieve_type::traits::coneSequence>& cone = 
@@ -144,15 +136,15 @@
   // check materials
   const int size = numCells;
   int* materialIds = (size > 0) ? new int[size] : 0;
-  const ALE::Obj<patch_label_type>& labelMaterials = 
-    topology->getLabel(patch, "material-id");
+  const ALE::Obj<label_type>& labelMaterials = 
+    mesh->getLabel("material-id");
   const int idDefault = -999;
 
   i = 0;
-  for(Mesh::topology_type::label_sequence::iterator e_iter = cells->begin();
+  for(Mesh::label_sequence::iterator e_iter = cells->begin();
       e_iter != cells->end();
       ++e_iter)
-    materialIds[i++] = topology->getValue(labelMaterials, *e_iter, idDefault);
+    materialIds[i++] = mesh->getValue(labelMaterials, *e_iter, idDefault);
   
   for (int iCell=0; iCell < numCells; ++iCell)
     CPPUNIT_ASSERT_EQUAL(data.materialIds[iCell], materialIds[iCell]);

Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.hh	2007-03-30 21:05:09 UTC (rev 6475)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.hh	2007-03-30 23:49:16 UTC (rev 6476)
@@ -39,9 +39,8 @@
   // PUBLIC TYPEDEFS ////////////////////////////////////////////////////
 public :
 
-  typedef ALE::Mesh Mesh;
-  typedef ALE::Mesh::sieve_type sieve_type;
-  typedef ALE::Mesh::topology_type topology_type;
+  typedef ALE::Field::Mesh Mesh;
+  typedef Mesh::sieve_type sieve_type;
 
   // PUBLIC METHODS /////////////////////////////////////////////////////
 public :
@@ -52,14 +51,14 @@
    *
    * @returns PETSc mesh
    */
-  ALE::Obj<ALE::Mesh>* createMesh(const MeshData& data);
+  ALE::Obj<Mesh>* createMesh(const MeshData& data);
 
   /** Check values in mesh against data.
    *
    * @param mesh PETSc mesh
    * @param data Mesh data
    */
-  void checkVals(const ALE::Obj<ALE::Mesh>& mesh,
+  void checkVals(const ALE::Obj<Mesh>& mesh,
 		 const MeshData& data);
 
 }; // class TestMeshIO



More information about the cig-commits mailing list