[cig-commits] r7688 - in short/3D/PyLith/trunk: . examples/3d/tet4 libsrc/faults libsrc/feassemble modulesrc/faults modulesrc/feassemble playpen/memcheck pylith pylith/bc pylith/faults pylith/feassemble pylith/feassemble/quadrature pylith/materials pylith/meshio pylith/problems unittests/pytests/bc unittests/pytests/faults unittests/pytests/feassemble unittests/pytests/materials

brad at geodynamics.org brad at geodynamics.org
Tue Jul 17 17:13:26 PDT 2007


Author: brad
Date: 2007-07-17 17:13:24 -0700 (Tue, 17 Jul 2007)
New Revision: 7688

Modified:
   short/3D/PyLith/trunk/TODO
   short/3D/PyLith/trunk/examples/3d/tet4/pylithapp.cfg
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
   short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
   short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
   short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
   short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
   short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
   short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
   short/3D/PyLith/trunk/playpen/memcheck/runtet4.cc
   short/3D/PyLith/trunk/pylith/PyLithApp.py
   short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py
   short/3D/PyLith/trunk/pylith/bc/Dirichlet.py
   short/3D/PyLith/trunk/pylith/faults/BruneSlipFn.py
   short/3D/PyLith/trunk/pylith/faults/EqKinSrc.py
   short/3D/PyLith/trunk/pylith/faults/Fault.py
   short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
   short/3D/PyLith/trunk/pylith/faults/SlipTimeFn.py
   short/3D/PyLith/trunk/pylith/feassemble/Constraint.py
   short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
   short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py
   short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature.py
   short/3D/PyLith/trunk/pylith/materials/Material.py
   short/3D/PyLith/trunk/pylith/meshio/SolutionIO.py
   short/3D/PyLith/trunk/pylith/problems/Explicit.py
   short/3D/PyLith/trunk/pylith/problems/Formulation.py
   short/3D/PyLith/trunk/pylith/problems/Implicit.py
   short/3D/PyLith/trunk/pylith/problems/Problem.py
   short/3D/PyLith/trunk/pylith/problems/TimeDependent.py
   short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichlet.py
   short/3D/PyLith/trunk/unittests/pytests/faults/TestBruneSlipFn.py
   short/3D/PyLith/trunk/unittests/pytests/faults/TestEqKinSrc.py
   short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py
   short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py
   short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py
   short/3D/PyLith/trunk/unittests/pytests/feassemble/TestQuadrature.py
   short/3D/PyLith/trunk/unittests/pytests/materials/TestMaterial.py
Log:
Split Python initialize() into preinitialize(), verifyConfiguration(), and initialize() to permit better error checking before simulation goes very far. Added tests for compatibility of quadrature scheme and cells for solid and cohesive cells.

Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/TODO	2007-07-18 00:13:24 UTC (rev 7688)
@@ -2,9 +2,13 @@
 CURRENT ISSUES
 ======================================================================
 
-  1. Need check to make sure quadrature scheme is compatible with
-  cells in mesh.
+  0. add shape to Quadrature
 
+  1. Memory analyzing.
+
+     Matt is working on buildTopology(), stratify(), restrict(), and
+     setValue().
+
   2. Add check to make sure every material in mesh has a material model.
 
   3. Need better error trapping when using LineParser. State of
@@ -16,7 +20,7 @@
 
   5. Add dependency diagram to manual.
 
-  6. Better default PETSc settings. (pc_type=asm?)
+  6. Switch to better default PETSc settings in examples. (pc_type=asm?)
 
   7. Experiment with different preconditioners and tabulate results
 

Modified: short/3D/PyLith/trunk/examples/3d/tet4/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/tet4/pylithapp.cfg	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/examples/3d/tet4/pylithapp.cfg	2007-07-18 00:13:24 UTC (rev 7688)
@@ -27,8 +27,8 @@
 
 [pylithapp.mesh_generator.importer]
 # Set filenames of mesh and pset files to import.
-filename_gmv = tet4_1000m_ascii.gmv
-filename_pset = tet4_1000m_ascii.pset
+filename_gmv = tet4_1000m_binary.gmv
+filename_pset = tet4_1000m_binary.pset
 
 # If using provided mesh or importing the mesh on a machine with a
 # different endian type than the one which created the mesh uncomment

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc	2007-07-18 00:13:24 UTC (rev 7688)
@@ -71,4 +71,45 @@
 } // integrateJacobian
   
 
+// ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+pylith::faults::FaultCohesiveDyn::verifyConfiguration(
+					      const ALE::Obj<Mesh>& mesh)
+{ // verifyConfiguration
+  assert(0 != _quadrature);
+
+  // check compatibility of mesh and quadrature scheme
+  const int dimension = mesh->getDimension()-1;
+  if (_quadrature->cellDim() != dimension) {
+    std::ostringstream msg;
+    msg << "Dimension of reference cell in quadrature scheme (" 
+	<< _quadrature->cellDim() 
+	<< ") does not match dimension of cells in mesh (" 
+	<< dimension << ") for fault '" << label()
+	<< "'.";
+    throw std::runtime_error(msg.str());
+  } // if
+  const int numCorners = _quadrature->numBasis();
+  const ALE::Obj<ALE::Mesh::label_sequence>& cells = 
+    mesh->getLabelStratum("material-id", id());
+  assert(!cells.isNull());
+  const Mesh::label_sequence::iterator cellsBegin = cells->begin();
+  const Mesh::label_sequence::iterator cellsEnd = cells->end();
+  const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
+  assert(!sieve.isNull());
+  for (Mesh::label_sequence::iterator c_iter=cellsBegin;
+       c_iter != cellsEnd;
+       ++c_iter) {
+    const int cellNumCorners = sieve->nCone(*c_iter, mesh->depth())->size();
+    if (3*numCorners != cellNumCorners) {
+      std::ostringstream msg;
+      msg << "Number of vertices in reference cell (" << numCorners 
+	  << ") is not compatible with number of vertices (" << cellNumCorners
+	  << ") in cohesive cell " << *c_iter << " for fault '"
+	  << label() << "'.";
+      throw std::runtime_error(msg.str());
+    } // if
+  } // for
+} // verifyConfiguration
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh	2007-07-18 00:13:24 UTC (rev 7688)
@@ -95,6 +95,12 @@
 			 topology::FieldsManager* const fields,
 			 const ALE::Obj<Mesh>& mesh);
   
+  /** Verify configuration is acceptable.
+   *
+   * @param mesh Finite-element mesh
+   */
+  void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+
   // PROTECTED METHODS //////////////////////////////////////////////////
 protected :
 

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2007-07-18 00:13:24 UTC (rev 7688)
@@ -619,5 +619,47 @@
   _needNewJacobian = false;
 } // integrateJacobian
   
+// ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+pylith::faults::FaultCohesiveKin::verifyConfiguration(
+					      const ALE::Obj<Mesh>& mesh)
+{ // verifyConfiguration
+  assert(0 != _quadrature);
 
+  // check compatibility of mesh and quadrature scheme
+  const int dimension = mesh->getDimension()-1;
+  if (_quadrature->cellDim() != dimension) {
+    std::ostringstream msg;
+    msg << "Dimension of reference cell in quadrature scheme (" 
+	<< _quadrature->cellDim() 
+	<< ") does not match dimension of cells in mesh (" 
+	<< dimension << ") for fault '" << label()
+	<< "'.";
+    throw std::runtime_error(msg.str());
+  } // if
+  const int numCorners = _quadrature->numBasis();
+  const ALE::Obj<ALE::Mesh::label_sequence>& cells = 
+    mesh->getLabelStratum("material-id", id());
+  assert(!cells.isNull());
+  const Mesh::label_sequence::iterator cellsBegin = cells->begin();
+  const Mesh::label_sequence::iterator cellsEnd = cells->end();
+  const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
+  assert(!sieve.isNull());
+  for (Mesh::label_sequence::iterator c_iter=cellsBegin;
+       c_iter != cellsEnd;
+       ++c_iter) {
+    const int cellNumCorners = sieve->nCone(*c_iter, mesh->depth())->size();
+    if (3*numCorners != cellNumCorners) {
+      std::ostringstream msg;
+      msg << "Number of vertices in reference cell (" << numCorners 
+	  << ") is not compatible with number of vertices (" << cellNumCorners
+	  << ") in cohesive cell " << *c_iter << " for fault '"
+	  << label() << "'.";
+      throw std::runtime_error(msg.str());
+    } // if
+  } // for
+} // verifyConfiguration
+
+
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh	2007-07-18 00:13:24 UTC (rev 7688)
@@ -118,6 +118,12 @@
 			 topology::FieldsManager* const fields,
 			 const ALE::Obj<Mesh>& mesh);
 
+  /** Verify configuration is acceptable.
+   *
+   * @param mesh Finite-element mesh
+   */
+  void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+
   // PROTECTED METHODS //////////////////////////////////////////////////
 protected :
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc	2007-07-18 00:13:24 UTC (rev 7688)
@@ -485,5 +485,60 @@
   } // for
 } // updateState
 
+// ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+pylith::feassemble::ElasticityExplicit::verifyConfiguration(
+						 const ALE::Obj<Mesh>& mesh)
+{ // verifyConfiguration
+  assert(0 != _quadrature);
+  assert(0 != _material);
 
+  const int dimension = mesh->getDimension();
+
+  // check compatibility of mesh and material
+  if (_material->dimension() != dimension) {
+    std::ostringstream msg;
+    msg << "Material '" << _material->label()
+	<< "' is incompatible with mesh.\n"
+	<< "Dimension of mesh: " << dimension
+	<< ", dimension of material: " << _material->dimension()
+	<< ".";
+    throw std::runtime_error(msg.str());
+  } // if
+
+  // check compatibility of mesh and quadrature scheme
+  if (_quadrature->cellDim() != dimension) {
+    std::ostringstream msg;
+    msg << "Quadrature is incompatible with cells for material '"
+	<< _material->label() << "'.\n"
+	<< "Dimension of mesh: " << dimension
+	<< ", dimension of quadrature: " << _quadrature->cellDim()
+	<< ".";
+    throw std::runtime_error(msg.str());
+  } // if
+  const int numCorners = _quadrature->numBasis();
+  const ALE::Obj<ALE::Mesh::label_sequence>& cells = 
+    mesh->getLabelStratum("material-id", _material->id());
+  assert(!cells.isNull());
+  const Mesh::label_sequence::iterator cellsEnd = cells->end();
+  const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
+  assert(!sieve.isNull());
+  for (Mesh::label_sequence::iterator c_iter=cells->begin();
+       c_iter != cellsEnd;
+       ++c_iter) {
+    const int cellNumCorners = sieve->nCone(*c_iter, mesh->depth())->size();
+    if (numCorners != cellNumCorners) {
+      std::ostringstream msg;
+      msg << "Quadrature is incompatible with cell in material '"
+	  << _material->label() << "'.\n"
+	  << "Cell " << *c_iter << " has " << cellNumCorners
+	  << " vertices but quadrature reference cell has "
+	  << numCorners << " vertices.";
+      throw std::runtime_error(msg.str());
+    } // if
+  } // for
+} // verifyConfiguration
+
+
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh	2007-07-18 00:13:24 UTC (rev 7688)
@@ -150,6 +150,12 @@
 		   const ALE::Obj<real_section_type>& field,
 		   const ALE::Obj<Mesh>& mesh);
 
+  /** Verify configuration is acceptable.
+   *
+   * @param mesh Finite-element mesh
+   */
+  void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+
 // PRIVATE METHODS //////////////////////////////////////////////////////
 private :
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc	2007-07-18 00:13:24 UTC (rev 7688)
@@ -369,6 +369,7 @@
 				   const ALE::Obj<real_section_type>& disp,
 				   const ALE::Obj<Mesh>& mesh)
 { // updateState
+  assert(0 != _quadrature);
   assert(0 != _material);
   assert(!disp.isNull());
 
@@ -441,6 +442,61 @@
 } // updateState
 
 // ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+pylith::feassemble::ElasticityImplicit::verifyConfiguration(
+						 const ALE::Obj<Mesh>& mesh)
+{ // verifyConfiguration
+  assert(0 != _quadrature);
+  assert(0 != _material);
+
+  const int dimension = mesh->getDimension();
+
+  // check compatibility of mesh and material
+  if (_material->dimension() != dimension) {
+    std::ostringstream msg;
+    msg << "Material '" << _material->label()
+	<< "' is incompatible with mesh.\n"
+	<< "Dimension of mesh: " << dimension
+	<< ", dimension of material: " << _material->dimension()
+	<< ".";
+    throw std::runtime_error(msg.str());
+  } // if
+
+  // check compatibility of mesh and quadrature scheme
+  if (_quadrature->cellDim() != dimension) {
+    std::ostringstream msg;
+    msg << "Quadrature is incompatible with cells for material '"
+	<< _material->label() << "'.\n"
+	<< "Dimension of mesh: " << dimension
+	<< ", dimension of quadrature: " << _quadrature->cellDim()
+	<< ".";
+    throw std::runtime_error(msg.str());
+  } // if
+  const int numCorners = _quadrature->numBasis();
+  const ALE::Obj<ALE::Mesh::label_sequence>& cells = 
+    mesh->getLabelStratum("material-id", _material->id());
+  assert(!cells.isNull());
+  const Mesh::label_sequence::iterator cellsEnd = cells->end();
+  const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
+  assert(!sieve.isNull());
+  for (Mesh::label_sequence::iterator c_iter=cells->begin();
+       c_iter != cellsEnd;
+       ++c_iter) {
+    const int cellNumCorners = sieve->nCone(*c_iter, mesh->depth())->size();
+    if (numCorners != cellNumCorners) {
+      std::ostringstream msg;
+      msg << "Quadrature is incompatible with cell in material '"
+	  << _material->label() << "'.\n"
+	  << "Cell " << *c_iter << " has " << cellNumCorners
+	  << " vertices but quadrature reference cell has "
+	  << numCorners << " vertices.";
+      throw std::runtime_error(msg.str());
+    } // if
+  } // for
+} // verifyConfiguration
+
+// ----------------------------------------------------------------------
 // Integrate elasticity term in residual for 1-D cells.
 void
 pylith::feassemble::ElasticityImplicit::_elasticityResidual1D(

Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh	2007-07-18 00:13:24 UTC (rev 7688)
@@ -150,6 +150,12 @@
 		   const ALE::Obj<real_section_type>& field,
 		   const ALE::Obj<Mesh>& mesh);
   
+  /** Verify configuration is acceptable.
+   *
+   * @param mesh Finite-element mesh
+   */
+  void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+
 // PRIVATE METHODS //////////////////////////////////////////////////////
 private :
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh	2007-07-18 00:13:24 UTC (rev 7688)
@@ -130,6 +130,13 @@
 		   const ALE::Obj<real_section_type>& field,
 		   const ALE::Obj<Mesh>& mesh);
 
+  /** Verify configuration is acceptable.
+   *
+   * @param mesh Finite-element mesh
+   */
+  virtual
+  void verifyConfiguration(const ALE::Obj<Mesh>& mesh) = 0;
+
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :
 

Modified: short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src	2007-07-18 00:13:24 UTC (rev 7688)
@@ -426,6 +426,37 @@
     return
 
 
+  def verifyConfiguration(self, mesh):
+    """
+    Verify compatibility of configuration settings.
+    """
+    # create shim for method 'verifyConfiguration'
+    #embed{ void FaultCohesiveKin_verifyConfiguration(void* objVptr, void* meshVptr)
+    try {
+      assert(0 != objVptr);
+      assert(0 != meshVptr);
+      ALE::Obj<ALE::Mesh>* mesh =
+        (ALE::Obj<ALE::Mesh>*) meshVptr;
+      ((pylith::faults::FaultCohesiveKin*) objVptr)->verifyConfiguration(*mesh);
+    } catch (const std::exception& err) {
+      PyErr_SetString(PyExc_RuntimeError,
+                      const_cast<char*>(err.what()));
+    } catch (const ALE::Exception& err) {
+      PyErr_SetString(PyExc_RuntimeError,
+                      const_cast<char*>(err.msg().c_str()));
+    } catch (...) {
+      PyErr_SetString(PyExc_RuntimeError,
+                      "Caught unknown C++ exception.");
+    } // try/catch
+    #}embed
+    if mesh.name != "pylith_topology_Mesh":
+      raise TypeError, \
+            "Argument 'mesh' must be extension module type 'Mesh'."
+    FaultCohesiveKin_verifyConfiguration(self.thisptr,
+                                   ptrFromHandle(mesh))
+    return
+
+
   property quadrature:
     def __set__(self, q):
       """

Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src	2007-07-18 00:13:24 UTC (rev 7688)
@@ -1094,6 +1094,37 @@
     return
 
 
+  def verifyConfiguration(self, mesh):
+    """
+    Verify compatibility of configuration settings.
+    """
+    # create shim for method 'verifyConfiguration'
+    #embed{ void Integrator_verifyConfiguration(void* objVptr, void* meshVptr)
+    try {
+      assert(0 != objVptr);
+      assert(0 != meshVptr);
+      ALE::Obj<ALE::Mesh>* mesh =
+        (ALE::Obj<ALE::Mesh>*) meshVptr;
+      ((pylith::feassemble::Integrator*) objVptr)->verifyConfiguration(*mesh);
+    } catch (const std::exception& err) {
+      PyErr_SetString(PyExc_RuntimeError,
+                      const_cast<char*>(err.what()));
+    } catch (const ALE::Exception& err) {
+      PyErr_SetString(PyExc_RuntimeError,
+                      const_cast<char*>(err.msg().c_str()));
+    } catch (...) {
+      PyErr_SetString(PyExc_RuntimeError,
+                      "Caught unknown C++ exception.");
+    } // try/catch
+    #}embed
+    if mesh.name != "pylith_topology_Mesh":
+      raise TypeError, \
+            "Argument 'mesh' must be extension module type 'Mesh'."
+    Integrator_verifyConfiguration(self.thisptr,
+                                   ptrFromHandle(mesh))
+    return
+
+
   def _createHandle(self):
     """
     Wrap pointer to C++ object in PyCObject.

Modified: short/3D/PyLith/trunk/playpen/memcheck/runtet4.cc
===================================================================
--- short/3D/PyLith/trunk/playpen/memcheck/runtet4.cc	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/playpen/memcheck/runtet4.cc	2007-07-18 00:13:24 UTC (rev 7688)
@@ -49,10 +49,10 @@
     iohandler.filenamePset(filenamePset);
     ALE::Obj<pylith::Mesh> mesh;
     iohandler.read(&mesh);
-
     spatialdata::geocoords::CSCart cs;
     cs.initialize();
 
+
     pylith::feassemble::Quadrature3D quadrature;
     const int cellDim = 3;
     const int numCorners = 4;

Modified: short/3D/PyLith/trunk/pylith/PyLithApp.py
===================================================================
--- short/3D/PyLith/trunk/pylith/PyLithApp.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/PyLithApp.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -81,8 +81,12 @@
       interfaces = self.problem.interfaces.bin
     mesh = self.mesher.create(interfaces)
 
-    # Initialize problem and then run
-    self.problem.initialize(mesh)
+    # Setup problem, verify configuration, and then initialize
+    self.problem.preinitialize(mesh)
+    self.problem.verifyConfiguration()
+    self.problem.initialize()
+
+    # Run problem and cleanup
     self.problem.run(self)
     self.problem.finalize()
     

Modified: short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -77,21 +77,36 @@
     return
 
 
-  def initialize(self, mesh):
+  def preinitialize(self, mesh):
     """
-    Initialize boundary condition.
+    Setup boundary condition.
     """
     self._createCppHandle()
-    
-    self.db.initialize()
     self.cppHandle.id = self.id
     self.cppHandle.label = self.label
-    self.cppHandle.db = self.db.cppHandle    
     self.mesh = mesh
-    self.cppHandle.initialize(mesh.cppHandle, mesh.coordsys.cppHandle)
     return
 
 
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    # :TODO:  Make sure mesh contains group with label.
+    return
+
+
+  def initialize(self):
+    """
+    Initialize boundary condition.
+    """    
+    self.db.initialize()
+    self.cppHandle.db = self.db.cppHandle    
+    self.cppHandle.initialize(self.mesh.cppHandle,
+                              self.mesh.coordsys.cppHandle)
+    return
+
+
   # PRIVATE METHODS ////////////////////////////////////////////////////
 
   def _configure(self):

Modified: short/3D/PyLith/trunk/pylith/bc/Dirichlet.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/Dirichlet.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/bc/Dirichlet.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -81,15 +81,21 @@
     return
 
 
-  def initialize(self, mesh):
+  def preinitialize(self, mesh):
     """
-    Initialize Dirichlet boundary condition.
+    Do pre-initialization setup.
     """
-    self._createCppHandle()
-
+    BoundaryCondition.preinitialize(self, mesh)
     self.cppHandle.fixedDOF = self.fixedDOF    
-    BoundaryCondition.initialize(self, mesh)
     return
+
+
+  def initialize(self):
+    """
+    Initialize Dirichlet boundary condition.
+    """
+    BoundaryCondition.initialize(self)
+    return
   
 
   # PRIVATE METHODS ////////////////////////////////////////////////////

Modified: short/3D/PyLith/trunk/pylith/faults/BruneSlipFn.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/BruneSlipFn.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/faults/BruneSlipFn.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -75,14 +75,19 @@
     return
 
 
+  def preinitialize(self):
+    """
+    Do pre-initialization setup.
+    """
+    self._createCppHandle()      
+    SlipTimeFn.preinitialize(self)
+    return
+
+
   def initialize(self):
     """
     Initialize.
     """
-    self._createCppHandle()
-      
-    SlipTimeFn.initialize(self)
-
     self.slip.initialize()
     self.slipTime.initialize()
     self.slipRate.initialize()

Modified: short/3D/PyLith/trunk/pylith/faults/EqKinSrc.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/EqKinSrc.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/faults/EqKinSrc.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -66,13 +66,29 @@
     return
 
 
+  def preinitialize(self):
+    """
+    Do pre-initialization setup.
+    """
+    self._createCppHandle()
+    self.slipfn.preinitialize()
+    self.cppHandle.slipfn = self.slipfn.cppHandle
+    return
+  
+
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    self.slipfn.verifyConfiguration()
+    return
+
+
   def initialize(self):
     """
     Initialize.
     """
-    self._createCppHandle()
     self.slipfn.initialize()
-    self.cppHandle.slipfn = self.slipfn.cppHandle
     return
 
 

Modified: short/3D/PyLith/trunk/pylith/faults/Fault.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/Fault.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/faults/Fault.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -82,7 +82,7 @@
                      "and materials)."
 
     label = pyre.inventory.str("label", default="")
-    label.meta['tip'] = "Name of material."
+    label.meta['tip'] = "Name of fault."
 
     upDir = pyre.inventory.list("up_dir", default=[0, 0, 1],
                                 validator=validateDir)
@@ -116,6 +116,7 @@
     """
     Component.__init__(self, name, facility="fault")
     self.cppHandle = None
+    self.mesh = None
     return
 
 
@@ -132,27 +133,46 @@
     return
 
 
-  def initialize(self, mesh):
+  def preinitialize(self, mesh):
     """
-    Initialize fault.
+    Setup fault.
     """
     self._createCppHandle()
+    self.cppHandle.id = self.id
+    self.cppHandle.label = self.label
+
+    self.mesh = mesh
     
-    self.quadrature.initialize()
-    self.matDB.initialize()
+    self.quadrature.preinitialize()    
+    self.cppHandle.quadrature = self.quadrature.cppHandle
+    return
+  
 
-    faultDim = mesh.dimension() - 1
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    faultDim = self.mesh.dimension() - 1
     if faultDim != self.quadrature.cell.cellDim:
       raise ValueError, \
             "Quadrature is incompatible with fault surface.\n" \
             "Dimensions for quadrature: %d, dimensions of fault: %d" % \
             (self.quadrature.cell.cellDim, faultDim)
 
+    # :TODO: Make sure mesh has group of vertices with label.
+    return
+  
+
+  def initialize(self):
+    """
+    Initialize fault.
+    """
+    self.quadrature.initialize()
+    self.matDB.initialize()
+
     assert(None != self.cppHandle)
-    self.cppHandle.id = self.id
-    self.cppHandle.label = self.label
-    self.cppHandle.quadrature = self.quadrature.cppHandle
-    self.cppHandle.initialize(mesh.cppHandle, mesh.coordsys.cppHandle,
+    self.cppHandle.initialize(self.mesh.cppHandle,
+                              self.mesh.coordsys.cppHandle,
                               self.upDir, self.normalDir,
                               self.matDB.cppHandle)
     return

Modified: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -63,18 +63,35 @@
     return
 
 
-  def initialize(self, mesh):
+  def preinitialize(self, mesh):
     """
+    Do pre-initialization setup.
+    """
+    self._info.log("Pre-initializing fault '%s'." % self.label)
+    FaultCohesive.preinitialize(self, mesh)
+    assert(None != self.cppHandle)
+    self.eqsrc.preinitialize()
+    self.cppHandle.eqsrc = self.eqsrc.cppHandle
+    return
+  
+
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    FaultCohesive.verifyConfiguration(self)
+    assert(None != self.cppHandle)
+    self.cppHandle.verifyConfiguration(self.mesh.cppHandle)
+    return
+
+
+  def initialize(self):
+    """
     Initialize cohesive elements.
     """
     self._info.log("Initializing fault '%s'." % self.label)
-    self._createCppHandle
-    
-    self.mesh = mesh
-    assert(None != self.cppHandle)
     self.eqsrc.initialize()
-    self.cppHandle.eqsrc = self.eqsrc.cppHandle
-    FaultCohesive.initialize(self, mesh)
+    FaultCohesive.initialize(self)
     return
 
 

Modified: short/3D/PyLith/trunk/pylith/faults/SlipTimeFn.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/SlipTimeFn.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/faults/SlipTimeFn.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -38,6 +38,20 @@
     return
 
 
+  def preinitialize(self):
+    """
+    Do pre-initialization setup.
+    """
+    return
+
+
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    return
+
+
   def initialize(self):
     """
     Initialize.

Modified: short/3D/PyLith/trunk/pylith/feassemble/Constraint.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/Constraint.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/feassemble/Constraint.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -48,10 +48,26 @@
     """
     Constructor.
     """
+    self.cppHandle = None
     self.mesh = None
     return
 
 
+  def preinitialize(self, mesh):
+    """
+    Setup constraint.
+    """
+    self.mesh = mesh
+    return
+
+
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    return
+
+
   def timeStep(self, dt):
     """
     Set time step for advancing from time t to time t+dt.

Modified: short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/Integrator.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/feassemble/Integrator.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -29,6 +29,7 @@
      not "integrateResidual" in attrs or \
      not "integrateJacobian" in attrs or \
      not "updateState" in attrs or \
+     not "verifyConfiguration" in attrs or \
      not "finalize" in attrs:
     result = False
   return result
@@ -54,25 +55,15 @@
     return
 
 
-  def setMesh(self, mesh):
+  def verifyConfiguration(self):
     """
-    Set mesh.
+    Verify compatibility of configuration.
     """
-    self.mesh = mesh
-    return
-  
-
-  def initQuadrature(self, quadrature):
-    """
-    Initialize quadrature.
-    """
     assert(None != self.cppHandle)
-    quadrature.initialize()
-    self.quadrature = quadrature
-    self.cppHandle.quadrature = self.quadrature.cppHandle
+    self.cppHandle.verifyConfiguration(self.mesh.cppHandle)
     return
-  
-  
+
+
   def timeStep(self, dt):
     """
     Set time step for advancing from time t to time t+dt.

Modified: short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -38,15 +38,32 @@
     return
 
 
-  def initMaterial(self, mesh, material):
+  def preinitialize(self, mesh, material):
     """
-    Initialize material properties.
+    Setup integrator.
     """
-    self._info.log("Initializing integrator material '%s'." % material.label)
-    material.initialize(mesh)
+    assert(None != self.cppHandle)
+
+    self.mesh = mesh
+    
+    material.preinitialize()
+
+    self.quadrature = material.quadrature
+    self.cppHandle.quadrature = self.quadrature.cppHandle
+
     self.material = material
     self.cppHandle.material = self.material.cppHandle
     return
   
+
+  def initialize(self):
+    """
+    Initialize material properties.
+    """
+    self._info.log("Initializing integrator for material '%s'." % \
+                   self.material.label)
+    self.material.initialize(self.mesh)
+    return
   
+  
 # End of file 

Modified: short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -67,9 +67,9 @@
     return
 
 
-  def initialize(self):
+  def preinitialize(self):
     """
-    Initialize C++ quadrature object.
+    Setup quadrature object.
     """
     self._createCppHandle()
     
@@ -95,6 +95,13 @@
     return
 
 
+  def initialize(self):
+    """
+    Initialize quadrature object.
+    """
+    return
+
+
   # PRIVATE METHODS ////////////////////////////////////////////////////
 
   def _configure(self):

Modified: short/3D/PyLith/trunk/pylith/materials/Material.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/Material.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/materials/Material.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -86,27 +86,41 @@
     return
 
 
+  def preinitialize(self):
+    """
+    Do pre-initialization setup.
+    """
+    self._createCppHandle()
+    self.cppHandle.id = self.id
+    self.cppHandle.label = self.label
+    self.quadrature.preinitialize()
+    return
+
+
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    if self.quadrature.spaceDim != self.dimension:
+        raise ValueError, \
+              "Quadrature scheme and material are incompatible.\n" \
+              "Dimension for quadrature: %d\n" \
+              "Dimension for material '%s': %d" % \
+              (self.quadrature.spaceDim, self.label, self.dimension)
+    
+    # :TODO: Make sure mesh contains material (need to account for the
+    # fact that any given processor may only have a subset of the
+    # materials)
+    return
+  
+
   def initialize(self, mesh):
     """
     Initialize material property manager.
     """
     self._info.log("Initializing material '%s'." % self.label)
-    self._createCppHandle()
-
-    if self.dimension != self.quadrature.cell.cellDim:
-      raise ValueError, \
-            "Quadrature is incompatible with material.\n" \
-            "Dimensions for quadrature: %d, dimensions for material: %d" % \
-            (self.quadrature.cell.cellDim, self.dimension)
-    if self.dimension != mesh.dimension():
-      raise ValueError, \
-            "Material is incompatible with mesh.\n" \
-            "Dimensions for mesh: %d, dimensions for material: %d" % \
-            (mesh.dimension(), self.dimension)
-
+    assert(None != self.cppHandle)
     self.db.initialize()
-    self.cppHandle.id = self.id
-    self.cppHandle.label = self.label
     self.cppHandle.db = self.db.cppHandle
     self.cppHandle.initialize(mesh.cppHandle, mesh.coordsys.cppHandle,
                               self.quadrature.cppHandle)

Modified: short/3D/PyLith/trunk/pylith/meshio/SolutionIO.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/SolutionIO.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/meshio/SolutionIO.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -80,6 +80,13 @@
     return
 
 
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    return
+
+
   def open(self, mesh):
     """
     Open files for solution.

Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -58,15 +58,13 @@
     return ElasticityExplicit()
 
 
-  def initialize(self, mesh, materials, boundaryConditions,
-                 interfaceConditions, dimension, dt):
+  def initialize(self, dimension, dt):
     """
     Initialize problem for explicit time integration.
     """
     from pyre.units.time import second
     t = 0.0*second
-    Formulation.initialize(self, mesh, materials, boundaryConditions,
-                           interfaceConditions, dimension, dt)
+    Formulation.initialize(self, dimension, dt)
 
     self._info.log("Creating other fields and matrices.")
     self.fields.addReal("dispTpdt")

Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -56,6 +56,7 @@
     from pylith.meshio.SingleOutput import SingleOutput
     output = pyre.inventory.facility("output", family="object_bin",
                                      factory=SingleOutput)
+
   
   # PUBLIC METHODS /////////////////////////////////////////////////////
 
@@ -72,42 +73,33 @@
     return
 
 
-  def initialize(self, mesh, materials, boundaryConditions,
-                 interfaceConditions, dimension, dt):
+  def preinitialize(self, mesh, materials, boundaryConditions,
+                    interfaceConditions):
     """
-    Create integrators for each element family.
+    Create integrator for each element family.
     """
     from pylith.feassemble.Integrator import implementsIntegrator
     from pylith.feassemble.Constraint import implementsConstraint
 
-    from pylith.topology.FieldsManager import FieldsManager
-    self.fields = FieldsManager(mesh)
-    self.mesh   = mesh
+    self.mesh = mesh
     self.integrators = []
     self.constraints = []
 
-    self._info.log("Initializing materials.")
+    self._info.log("Pre-initializing materials.")
     for material in materials.bin:
-      if material.quadrature.spaceDim != dimension:
-        raise ValueError, \
-              "Spatial dimension of problem is '%d' but quadrature " \
-              "for material '%s' is for spatial dimension '%d'." % \
-              (dimension, material.label, material.quadrature.spaceDim)
       integrator = self.elasticityIntegrator()
       if not implementsIntegrator(integrator):
         raise TypeError, \
               "Could not use '%s' as an integrator for material '%s'. " \
               "Functionality missing." % (integrator.name, material.label)
-      integrator.setMesh(mesh)
-      integrator.initQuadrature(material.quadrature)
-      integrator.initMaterial(mesh, material)
+      integrator.preinitialize(mesh, material)
       self.integrators.append(integrator)
       self._info.log("Added elasticity integrator for material '%s'." % \
                      material.label)
 
-    self._info.log("Initializing boundary conditions.")
+    self._info.log("Pre-initializing boundary conditions.")
     for bc in boundaryConditions.bin:
-      bc.initialize(mesh)
+      bc.preinitialize(mesh)
       foundType = False
       if implementsIntegrator(bc):
         foundType = True
@@ -124,9 +116,9 @@
               "Could not determine whether boundary condition '%s' is an " \
               "integrator or a constraint." % bc.name
 
-    self._info.log("Initializing interior interfaces.")
+    self._info.log("Pre-initializing interior interfaces.")
     for ic in interfaceConditions.bin:
-      ic.initialize(mesh)
+      ic.preinitialize(mesh)
       foundType = False
       if implementsIntegrator(ic):
         foundType = True
@@ -142,10 +134,40 @@
         raise TypeError, \
               "Could not determine whether interface condition '%s' is an " \
               "integrator or a constraint." % ic.name
+    return
 
+
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    for integrator in self.integrators:
+      integrator.verifyConfiguration()
+    for constraint in self.constraints:
+      constraint.verifyConfiguration()
+    for output in self.output.bin:
+      output.verifyConfiguration()
+    return
+  
+
+  def initialize(self, dimension, dt):
+    """
+    Create integrators for each element family.
+    """
+    from pylith.topology.FieldsManager import FieldsManager
+    self.fields = FieldsManager(self.mesh)
+
+    self._info.log("Initializing integrators.")
+    for integrator in self.integrators:
+      integrator.initialize()
+
+    self._info.log("Initializing constraints.")
+    for constraint in self.constraints:
+      constraint.initialize()
+
     self._info.log("Setting up solution output.")
     for output in self.output.bin:
-      output.open(mesh)
+      output.open(self.mesh)
       output.writeTopology()
 
     self._info.log("Creating solution field.")

Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -81,20 +81,18 @@
     return ElasticityImplicit()
 
 
-  def initialize(self, mesh, materials, boundaryConditions,
-                 interfaceConditions, dimension, dt):
+  def initialize(self, dimension, dt):
     """
     Initialize problem for implicit time integration.
     """
-    Formulation.initialize(self, mesh, materials, boundaryConditions,
-                           interfaceConditions, dimension, dt)
+    Formulation.initialize(self, dimension, dt)
 
     self._info.log("Creating other fields and matrices.")
     self.fields.addReal("dispIncr")
     self.fields.addReal("residual")
     self.fields.copyLayout("dispTBctpdt")
-    self.jacobian = mesh.createMatrix(self.fields.getSolution())
-    self.solver.initialize(mesh, self.fields.getSolution())
+    self.jacobian = self.mesh.createMatrix(self.fields.getSolution())
+    self.solver.initialize(self.mesh, self.fields.getSolution())
 
     # Initial time step solves for total displacement field, not increment
     for constraint in self.constraints:

Modified: short/3D/PyLith/trunk/pylith/problems/Problem.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Problem.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/problems/Problem.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -73,11 +73,33 @@
     return
 
 
-  def initialize(self, mesh):
+  def preinitialize(self, mesh):
     """
     Setup integrators for each element family (material/quadrature,
     bc/quadrature, etc.).
     """
+    raise NotImplementedError, "preinitialize() not implemented."
+    return
+
+
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    for material in self.materials.bin:
+      material.verifyConfiguration()
+    for bc in self.bc.bin:
+      bc.verifyConfiguration()
+    for interface in self.interfaces.bin:
+      interface.verifyConfiguration()
+    return
+  
+
+  def initialize(self):
+    """
+    Initialize integrators for each element family (material/quadrature,
+    bc/quadrature, etc.).
+    """
     raise NotImplementedError, "initialize() not implemented."
     return
 

Modified: short/3D/PyLith/trunk/pylith/problems/TimeDependent.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/TimeDependent.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/pylith/problems/TimeDependent.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -84,24 +84,49 @@
     return
 
 
-  def initialize(self, mesh):
+  def preinitialize(self, mesh):
     """
     Setup integrators for each element family (material/quadrature,
     bc/quadrature, etc.).
     """
-    self._info.log("Initializing problem.")
+    self._info.log("Pre-initializing problem.")
+    self.mesh = mesh
+    self.formulation.preinitialize(mesh, self.materials, self.bc,
+                                   self.interfaces)
+    return
 
-    if self.dimension != mesh.dimension():
+
+  def verifyConfiguration(self):
+    """
+    Verify compatibility of configuration.
+    """
+    self._info.log("Verifying compatibility of problem configuration.")
+    if self.dimension != self.mesh.dimension():
+      raise ValueError, \
+            "Spatial dimension of problem is '%d' but mesh contains cells " \
+            "for spatial dimension '%d'." % \
+            (self.dimension, mesh.dimension)
+    for material in self.materials.bin:
+      if material.quadrature.spaceDim != self.dimension:
         raise ValueError, \
-              "Spatial dimension of problem is '%d' but mesh contains cells " \
-              "for spatial dimension '%d'." % \
-              (self.dimension, mesh.dimension)
-    self.mesh = mesh
-    self.formulation.initialize(mesh, self.materials, self.bc,
-                                self.interfaces, self.dimension, self.dt)
+              "Spatial dimension of problem is '%d' but quadrature " \
+              "for material '%s' is for spatial dimension '%d'." % \
+              (self.dimension, material.label, material.quadrature.spaceDim)
+    Problem.verifyConfiguration(self)
+    self.formulation.verifyConfiguration()
     return
+  
 
+  def initialize(self):
+    """
+    Setup integrators for each element family (material/quadrature,
+    bc/quadrature, etc.).
+    """
+    self._info.log("Initializing problem.")
+    self.formulation.initialize(self.dimension, self.dt)
+    return
 
+
   def run(self, app):
     """
     Solve time dependent problem.

Modified: short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichlet.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichlet.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichlet.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -191,7 +191,8 @@
     importer.coordsys = cs
     mesh = importer.read(debug=False, interpolate=False)
     
-    bc.initialize(mesh)
+    bc.preinitialize(mesh)
+    bc.initialize()
     return (mesh, bc)
 
 

Modified: short/3D/PyLith/trunk/unittests/pytests/faults/TestBruneSlipFn.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/faults/TestBruneSlipFn.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/unittests/pytests/faults/TestBruneSlipFn.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -63,6 +63,7 @@
     slipFn.slip = dbFinalSlip
     slipFn.slipTime = dbSlipTime
     slipFn.slipRate = dbPeakRate
+    slipFn.preinitialize()
     slipFn.initialize()
     return
 

Modified: short/3D/PyLith/trunk/unittests/pytests/faults/TestEqKinSrc.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/faults/TestEqKinSrc.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/unittests/pytests/faults/TestEqKinSrc.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -67,6 +67,7 @@
 
     eqsrc = EqKinSrc()
     eqsrc.slipfn = slipfn
+    eqsrc.preinitialize()
     eqsrc.initialize()
     return
 

Modified: short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -282,7 +282,8 @@
     fault.matDB = dbMat
     fault.timeStep(dt)
     fault.adjustTopology(mesh)
-    fault.initialize(mesh)
+    fault.preinitialize(mesh)
+    fault.initialize()
 
     # Setup fields
     from pylith.topology.FieldsManager import FieldsManager

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -36,45 +36,53 @@
     return
     
 
-  def test_initQuadrature(self):
+  def test_preinitialize(self):
     """
-    Test initQuadrature().
+    Test preiniitlaize().
     """
-    from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
-    q = Quadrature2D()
-    minJacobian = 4.0e-02;
-    q.minJacobian = minJacobian
-    from pylith.feassemble.FIATSimplex import FIATSimplex
-    q.cell = FIATSimplex()
-    q.cell.shape = "triangle"
-    q.cell.order = 1
-    q.cell.degree = 1
-
-    integrator = ElasticityExplicit()
-    integrator.initQuadrature(q)
-    self.assertEqual(minJacobian, integrator.quadrature.minJacobian)
-    return
-    
-
-  def test_setMesh(self):
-    """
-    Test setMesh().
-    """
+    # Setup mesh
     cs = CSCart()
     cs.spaceDim = 2
-    
     from pylith.meshio.MeshIOAscii import MeshIOAscii
     importer = MeshIOAscii()
     importer.filename = "data/tri3.mesh"
     importer.coordsys = cs
     mesh = importer.read(debug=False, interpolate=False)
 
+    # Setup material
+    from pylith.feassemble.FIATSimplex import FIATSimplex
+    cell = FIATSimplex()
+    cell.shape = "triangle"
+    cell.degree = 1
+    cell.order = 1
+    from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
+    quadrature = Quadrature2D()
+    quadrature.cell = cell
+    minJacobian = 4.0e-02;
+    quadrature.minJacobian = minJacobian
+    
+    from spatialdata.spatialdb.SimpleDB import SimpleDB
+    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+    iohandler = SimpleIOAscii()
+    iohandler.filename = "data/elasticplanestrain.spatialdb"
+    db = SimpleDB()
+    db.label = "elastic plane strain"
+    db.iohandler = iohandler
+
+    from pylith.materials.ElasticPlaneStrain import ElasticPlaneStrain
+    material = ElasticPlaneStrain()
+    material.id = 0
+    material.label = "elastic plane strain"
+    material.db = db
+    material.quadrature = quadrature
+
     integrator = ElasticityExplicit()
-    integrator.setMesh(mesh)
+    integrator.preinitialize(mesh, material)
     self.assertEqual(mesh, integrator.mesh)
+    self.assertEqual(minJacobian, integrator.quadrature.minJacobian)
     return
+    
 
-  
   def test_timeStep(self):
     """
     Test timeStep().
@@ -236,9 +244,8 @@
 
     # Setup integrator
     integrator = ElasticityExplicit()
-    integrator.setMesh(mesh)
-    integrator.initQuadrature(material.quadrature)
-    integrator.initMaterial(mesh, material)
+    integrator.preinitialize(mesh, material)
+    integrator.initialize()
     integrator.timeStep(dt)
 
     # Setup fields

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -36,45 +36,53 @@
     return
     
 
-  def test_initQuadrature(self):
+  def test_preinitialize(self):
     """
-    Test initQuadrature().
+    Test preiniitlaize().
     """
-    from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
-    q = Quadrature2D()
-    minJacobian = 4.0e-02;
-    q.minJacobian = minJacobian
-    from pylith.feassemble.FIATSimplex import FIATSimplex
-    q.cell = FIATSimplex()
-    q.cell.shape = "triangle"
-    q.cell.order = 1
-    q.cell.degree = 1
-
-    integrator = ElasticityImplicit()
-    integrator.initQuadrature(q)
-    self.assertEqual(minJacobian, integrator.quadrature.minJacobian)
-    return
-    
-
-  def test_setMesh(self):
-    """
-    Test setMesh().
-    """
+    # Setup mesh
     cs = CSCart()
     cs.spaceDim = 2
-    
     from pylith.meshio.MeshIOAscii import MeshIOAscii
     importer = MeshIOAscii()
     importer.filename = "data/tri3.mesh"
     importer.coordsys = cs
     mesh = importer.read(debug=False, interpolate=False)
 
+    # Setup material
+    from pylith.feassemble.FIATSimplex import FIATSimplex
+    cell = FIATSimplex()
+    cell.shape = "triangle"
+    cell.degree = 1
+    cell.order = 1
+    from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
+    quadrature = Quadrature2D()
+    quadrature.cell = cell
+    minJacobian = 4.0e-02;
+    quadrature.minJacobian = minJacobian
+    
+    from spatialdata.spatialdb.SimpleDB import SimpleDB
+    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+    iohandler = SimpleIOAscii()
+    iohandler.filename = "data/elasticplanestrain.spatialdb"
+    db = SimpleDB()
+    db.label = "elastic plane strain"
+    db.iohandler = iohandler
+
+    from pylith.materials.ElasticPlaneStrain import ElasticPlaneStrain
+    material = ElasticPlaneStrain()
+    material.id = 0
+    material.label = "elastic plane strain"
+    material.db = db
+    material.quadrature = quadrature
+
     integrator = ElasticityImplicit()
-    integrator.setMesh(mesh)
+    integrator.preinitialize(mesh, material)
     self.assertEqual(mesh, integrator.mesh)
+    self.assertEqual(minJacobian, integrator.quadrature.minJacobian)
     return
+    
 
-  
   def test_timeStep(self):
     """
     Test timeStep().
@@ -235,9 +243,8 @@
 
     # Setup integrator
     integrator = ElasticityImplicit()
-    integrator.setMesh(mesh)
-    integrator.initQuadrature(material.quadrature)
-    integrator.initMaterial(mesh, material)
+    integrator.preinitialize(mesh, material)
+    integrator.initialize()
     integrator.timeStep(dt)
 
     # Setup fields

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestQuadrature.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestQuadrature.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestQuadrature.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -93,6 +93,7 @@
     quadrature = Quadrature1D()
     quadrature.cell = cell
 
+    quadrature.preinitialize()
     quadrature.initialize()
 
     from pylith.utils.testarray import test_double

Modified: short/3D/PyLith/trunk/unittests/pytests/materials/TestMaterial.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/TestMaterial.py	2007-07-17 23:06:19 UTC (rev 7687)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/TestMaterial.py	2007-07-18 00:13:24 UTC (rev 7688)
@@ -41,7 +41,7 @@
     cell.degree = 1
     quadrature.cell = cell
     quadrature.minJacobian = 1.0e-4
-    quadrature.initialize()
+    quadrature.preinitialize()
     material.quadrature = quadrature
 
     from spatialdata.spatialdb.SimpleDB import SimpleDB
@@ -65,6 +65,7 @@
     importer.coordsys = cs
     mesh = importer.read(debug=False, interpolate=False)
     
+    material.preinitialize()
     material.initialize(mesh)
 
     # We should really add something here to check to make sure things



More information about the cig-commits mailing list