[cig-commits] r15059 - in short/3D/PyLith/trunk: . examples/3d/hex8 libsrc/faults modulesrc/faults pylith/faults tests/1d/line2 unittests/libtests/faults unittests/pytests/faults

brad at geodynamics.org brad at geodynamics.org
Tue May 26 14:59:32 PDT 2009


Author: brad
Date: 2009-05-26 14:59:31 -0700 (Tue, 26 May 2009)
New Revision: 15059

Modified:
   short/3D/PyLith/trunk/README
   short/3D/PyLith/trunk/examples/3d/hex8/dislocation.cfg
   short/3D/PyLith/trunk/examples/3d/hex8/savageprescott.cfg
   short/3D/PyLith/trunk/libsrc/faults/Fault.hh
   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/modulesrc/faults/Fault.i
   short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i
   short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveKin.i
   short/3D/PyLith/trunk/pylith/faults/Fault.py
   short/3D/PyLith/trunk/tests/1d/line2/dislocation.cfg
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
   short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py
Log:
Removed fault conditioning (no longer needed with nondimensionalization).

Modified: short/3D/PyLith/trunk/README
===================================================================
--- short/3D/PyLith/trunk/README	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/README	2009-05-26 21:59:31 UTC (rev 15059)
@@ -64,8 +64,6 @@
 user sets the component to a facility. This applies to gravity,
 initial stresses, initial strains, and initial state variables.
 
-STILL ON THE TODO LIST
-
 (7) Nondimensionalization of the problem eliminates the need to
 condition the fault constraints. The "mat_db" facility was removed.
 

Modified: short/3D/PyLith/trunk/examples/3d/hex8/dislocation.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/dislocation.cfg	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/examples/3d/hex8/dislocation.cfg	2009-05-26 21:59:31 UTC (rev 15059)
@@ -71,7 +71,6 @@
 label = 10
 quadrature.cell = pylith.feassemble.FIATLagrange
 quadrature.cell.dimension = 2
-mat_db.iohandler.filename = mat_elastic.spatialdb
 
 [pylithapp.timedependent.interfaces.fault.eq_srcs.rupture.slip_function]
 slip.iohandler.filename = finalslip.spatialdb

Modified: short/3D/PyLith/trunk/examples/3d/hex8/savageprescott.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/savageprescott.cfg	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/examples/3d/hex8/savageprescott.cfg	2009-05-26 21:59:31 UTC (rev 15059)
@@ -90,7 +90,6 @@
 label = 10
 quadrature.cell = pylith.feassemble.FIATLagrange
 quadrature.cell.dimension = 2
-mat_db.iohandler.filename = mat_elastic.spatialdb
 
 output.vertex_info_fields = [normal_dir,strike_dir,dip_dir,final_slip_creep,final_slip_one,final_slip_two,final_slip_three,slip_time_creep,slip_time_one,slip_time_two,slip_time_three]
 

Modified: short/3D/PyLith/trunk/libsrc/faults/Fault.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Fault.hh	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/libsrc/faults/Fault.hh	2009-05-26 21:59:31 UTC (rev 15059)
@@ -92,14 +92,11 @@
    * @param normalDir General preferred direction for fault normal
    *   (used to pick which of two possible normal directions for
    *   interface; only applies to fault surfaces in a 3-D domain).
-   * @param matDB Database of bulk elastic properties for fault region
-   *   (used to improve conditioning of Jacobian matrix)
    */
   virtual
   void initialize(const topology::Mesh& mesh,
 		  const double upDir[3],
-		  const double normalDir[3],
-		  spatialdata::spatialdb::SpatialDB* matDB) = 0;
+		  const double normalDir[3]) = 0;
 
   /** Get mesh associated with fault fields.
    *

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc	2009-05-26 21:59:31 UTC (rev 15059)
@@ -38,8 +38,7 @@
 void
 pylith::faults::FaultCohesiveDyn::initialize(const topology::Mesh& mesh,
 					     const double upDir[3],
-					     const double normalDir[3],
-					     spatialdata::spatialdb::SpatialDB* matDB)
+					     const double normalDir[3])
 { // initialize
   throw std::logic_error("FaultCohesiveDyn::initialize() not implemented.");
 } // initialize

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh	2009-05-26 21:59:31 UTC (rev 15059)
@@ -56,13 +56,10 @@
    * @param normalDir General preferred direction for fault normal
    *   (used to pick which of two possible normal directions for
    *   interface; only applies to fault surfaces in a 3-D domain).
-   * @param matDB Database of bulk elastic properties for fault region
-   *   (used to improve conditioning of Jacobian matrix)
    */
   void initialize(const topology::Mesh& mesh,
 		  const double upDir[3],
-		  const double normalDir[3],
-		  spatialdata::spatialdb::SpatialDB* matDB);
+		  const double normalDir[3]);
 
   /** Integrate contribution of cohesive cells to residual term.
    *

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2009-05-26 21:59:31 UTC (rev 15059)
@@ -40,7 +40,6 @@
 #include <stdexcept> // USES std::runtime_error
 
 //#define PRECOMPUTE_GEOMETRY
-//#define USE_CONDITIONING
 
 // ----------------------------------------------------------------------
 typedef pylith::topology::Mesh::SieveMesh SieveMesh;
@@ -86,8 +85,7 @@
 void
 pylith::faults::FaultCohesiveKin::initialize(const topology::Mesh& mesh,
 					     const double upDir[3],
-					     const double normalDir[3],
-					     spatialdata::spatialdb::SpatialDB* matDB)
+					     const double normalDir[3])
 { // initialize
   assert(0 != upDir);
   assert(0 != normalDir);
@@ -134,12 +132,6 @@
   cumSlip.vectorFieldType(topology::FieldBase::VECTOR);
   cumSlip.scale(_normalizer->lengthScale());
 
-#if defined(USE_CONDITIONING)
-  // Setup pseudo-stiffness of cohesive cells to improve conditioning
-  // of Jacobian matrix
-  _calcConditioning(cs, matDB);
-#endif
-
   const ALE::Obj<SieveSubMesh::label_sequence>& cells = 
     faultSieveMesh->heightStratum(0);
   assert(!cells.isNull());
@@ -192,9 +184,6 @@
 
   // Allocate vectors for cell values
   double_array orientationCell(numConstraintVert*orientationSize);
-#if defined(USE_CONDITIONING)
-  double_array stiffnessCell(numConstraintVert);
-#endif
   double_array dispTCell(numCorners*spaceDim);
   double_array residualCell(numCorners*spaceDim);
 
@@ -228,15 +217,6 @@
 						     orientationCell.size(),
 						     &orientationCell[0]);
 
-#if defined(USE_CONDITIONING)
-  const ALE::Obj<RealSection>& stiffnessSection = 
-    _fields->get("pseudostiffness").section();
-  assert(!stiffnessSection.isNull());
-  topology::Mesh::RestrictVisitor stiffnessVisitor(*stiffnessSection,
-						   stiffnessCell.size(),
-						   &stiffnessCell[0]);
-#endif
-
   const ALE::Obj<RealSection>& areaSection = 
     _fields->get("area").section();
   assert(!areaSection.isNull());
@@ -294,12 +274,6 @@
     orientationVisitor.clear();
     faultSieveMesh->restrictClosure(c_fault, orientationVisitor);
     
-#if defined(USE_CONDITIONING)
-    // Get pseudo stiffness at fault cell's vertices.
-    stiffnessVisitor.clear();
-    faultSieveMesh->restrictClosure(c_fault, stiffnessVisitor);
-#endif
-    
     // Get area at fault cell's vertices.
     areaVisitor.clear();
     faultSieveMesh->restrictClosure(c_fault, areaVisitor);
@@ -315,15 +289,8 @@
       const int indexJ = iConstraint +   numConstraintVert;
       const int indexK = iConstraint + 2*numConstraintVert;
 
-#if defined(USE_CONDITIONING)
-      const double pseudoStiffness = stiffnessCell[iConstraint];
       assert(areaCell[iConstraint] > 0);
-      const double wt = pseudoStiffness * 
-	areaVertexCell[iConstraint] / areaCell[iConstraint];
-#else
-      assert(areaCell[iConstraint] > 0);
       const double wt = areaVertexCell[iConstraint] / areaCell[iConstraint];
-#endif
       
       // Get orientation at constraint vertex
       const double* orientationVertex = 
@@ -490,9 +457,6 @@
   const int numCorners = 3*numConstraintVert; // cohesive cell
   double_array matrixCell(numCorners*spaceDim * numCorners*spaceDim);
   double_array orientationCell(numConstraintVert*orientationSize);
-#if defined(USE_CONDITIONING)
-  double_array stiffnessCell(numConstraintVert);
-#endif
 
   // Get section information
   const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
@@ -506,15 +470,6 @@
 						     orientationCell.size(),
 						     &orientationCell[0]);
 
-#if defined(USE_CONDITIONING)
-  const ALE::Obj<RealSection>& stiffnessSection = 
-    _fields->get("pseudostiffness").section();
-  assert(!stiffnessSection.isNull());
-  topology::Mesh::RestrictVisitor stiffnessVisitor(*stiffnessSection,
-						   stiffnessCell.size(),
-						   &stiffnessCell[0]);
-#endif
-
 #if 0 // DEBUGGING
   // Check that fault cells match cohesive cells
   ALE::ISieveVisitor::PointRetriever<sieve_type> cV(std::max(1, mesh->getSieve()->getMaxConeSize()));
@@ -566,12 +521,6 @@
     orientationVisitor.clear();
     faultSieveMesh->restrictClosure(c_fault, orientationVisitor);
 
-#if defined(USE_CONDITIONING)
-    // Get pseudo stiffness at fault cell's vertices.
-    stiffnessVisitor.clear();
-    faultSieveMesh->restrictClosure(c_fault, stiffnessVisitor);
-#endif
-    
     for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
       // Blocks in cell matrix associated with normal cohesive
       // vertices i and j and constraint vertex k
@@ -584,25 +533,13 @@
 	&orientationCell[iConstraint*orientationSize];
       assert(0 != orientationVertex);
 
-#if defined(USE_CONDITIONING)
-      const double stiffnessVertex = stiffnessCell[iConstraint];
-#endif
-
-      // Scale orientation information by pseudo-stiffness to bring
-      // constraint forces in solution vector to the same order of
-      // magnitude as the displacements to prevent ill-conditioning
-
       // Entries associated with constraint forces applied at node i
       for (int iDim=0; iDim < spaceDim; ++iDim)
 	for (int kDim=0; kDim < spaceDim; ++kDim) {
 	  const int row = indexI*spaceDim+iDim;
 	  const int col = indexK*spaceDim+kDim;
 	  matrixCell[row*numCorners*spaceDim+col] =
-#if defined(USE_CONDITIONING)
-	    -orientationVertex[kDim*spaceDim+iDim]*stiffnessVertex;
-#else
 	    -orientationVertex[kDim*spaceDim+iDim];
-#endif
 	  matrixCell[col*numCorners*spaceDim+row] =
 	    -orientationVertex[kDim*spaceDim+iDim];
 	} // for
@@ -613,11 +550,7 @@
 	  const int row = indexJ*spaceDim+jDim;
 	  const int col = indexK*spaceDim+kDim;
 	  matrixCell[row*numCorners*spaceDim+col] =
-#if defined(USE_CONDITIONING)
-	    orientationVertex[kDim*spaceDim+jDim]*stiffnessVertex;
-#else
 	    orientationVertex[kDim*spaceDim+jDim];
-#endif
 	  matrixCell[col*numCorners*spaceDim+row] =
 	    orientationVertex[kDim*spaceDim+jDim];
 	} // for
@@ -1006,84 +939,6 @@
 
 // ----------------------------------------------------------------------
 void
-pylith::faults::FaultCohesiveKin::_calcConditioning(
-				 const spatialdata::geocoords::CoordSys* cs,
-				 spatialdata::spatialdb::SpatialDB* matDB)
-{ // _calcConditioning
-  assert(0 != cs);
-  assert(0 != matDB);
-  assert(0 != _faultMesh);
-  assert(0 != _fields);
-
-  const int spaceDim = cs->spaceDim();
-
-  // Get vertices in fault mesh.
-  const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
-  assert(!faultSieveMesh.isNull());
-  const ALE::Obj<SieveSubMesh::label_sequence>& vertices = 
-    faultSieveMesh->depthStratum(0);
-  const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
-  const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
-  
-  // Allocate stiffness field.
-  _fields->add("pseudostiffness", "pseudostiffness");
-  topology::Field<topology::SubMesh>& stiffness = _fields->get("pseudostiffness");
-  stiffness.newSection(topology::FieldBase::VERTICES_FIELD, 1);
-  stiffness.allocate();
-  stiffness.zero();
-  const ALE::Obj<RealSection>& stiffnessSection = stiffness.section();
-  assert(!stiffnessSection.isNull());
-
-  // Setup queries of physical properties.
-  matDB->open();
-  const char* stiffnessVals[] = { "density", "vs" };
-  const int numStiffnessVals = 2;
-  matDB->queryVals(stiffnessVals, numStiffnessVals);
-  
-  // Get section containing coordinates of vertices
-  const ALE::Obj<RealSection>& coordinates = 
-    faultSieveMesh->getRealSection("coordinates");
-  assert(!coordinates.isNull());
-
-  // Get dimensional scales.
-  assert(0 != _normalizer);
-  const double pressureScale = _normalizer->pressureScale();
-  const double lengthScale = _normalizer->lengthScale();
-
-  double_array matprops(numStiffnessVals);
-  double_array coordsVertex(spaceDim);
-  int count = 0;
-  
-  // Set values in orientation section.
-  for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
-       v_iter != verticesEnd;
-       ++v_iter, ++count) {
-    coordinates->restrictPoint(*v_iter, &coordsVertex[0], coordsVertex.size());
-    _normalizer->dimensionalize(&coordsVertex[0], coordsVertex.size(), lengthScale);
-    int err = matDB->query(&matprops[0], numStiffnessVals, &coordsVertex[0], 
-			   coordsVertex.size(), cs);
-    if (err) {
-      std::ostringstream msg;
-      msg << "Could not find material properties at (";
-      for (int i=0; i < spaceDim; ++i)
-	msg << "  " << coordsVertex[i];
-      msg << ") using spatial database " << matDB->label() << ".";
-      throw std::runtime_error(msg.str());
-    } // if
-    
-    const double density = matprops[0];
-    const double vs = matprops[1];
-    const double mu = density * vs*vs;
-    const double muN = _normalizer->nondimensionalize(mu, pressureScale);
-    stiffnessSection->updatePoint(*v_iter, &muN);
-  } // for
-  PetscLogFlops(count * 2);
-
-  matDB->close();
-} // _calcConditioning
-
-// ----------------------------------------------------------------------
-void
 pylith::faults::FaultCohesiveKin::_calcArea(void)
 { // _calcArea
   assert(0 != _faultMesh);
@@ -1210,11 +1065,6 @@
   const SieveSubMesh::label_sequence::iterator fverticesEnd = fvertices->end();
 
   // Get sections.
-#if defined(USE_CONDITIONING)
-  const ALE::Obj<RealSection>& stiffnessSection = 
-    _fields->get("pseudostiffness").section();
-  assert(!stiffnessSection.isNull());
-#endif
   const ALE::Obj<RealSection>& areaSection = 
     _fields->get("area").section();
   assert(!areaSection.isNull());
@@ -1272,27 +1122,15 @@
       const int vertexFault = renumbering[*v_iter];
       assert(fiberDim == dispTSection->getFiberDimension(vertexMesh));
       assert(fiberDim == tractionsSection->getFiberDimension(vertexFault));
-#if defined(USE_CONDITIONING)
-      assert(1 == stiffnessSection->getFiberDimension(vertexFault));
-#endif
       assert(1 == areaSection->getFiberDimension(vertexFault));
 
       const double* dispTVertex = dispTSection->restrictPoint(vertexMesh);
       assert(0 != dispTVertex);
-#if defined(USE_CONDITIONING)
-      const double* stiffnessVertex = stiffnessSection->restrictPoint(vertexFault);
-      assert(0 != stiffnessVertex);
-#endif
       const double* areaVertex = areaSection->restrictPoint(vertexFault);
       assert(0 != areaVertex);
 
-#if defined(USE_CONDITIONING)
-      const double scale = stiffnessVertex[0] / areaVertex[0];
-#else
-      const double scale = 1.0 / areaVertex[0];
-#endif
       for (int i=0; i < fiberDim; ++i)
-	tractionsVertex[i] = dispTVertex[i] * scale;
+	tractionsVertex[i] = dispTVertex[i] / areaVertex[0];
 
       tractionsSection->updatePoint(vertexFault, &tractionsVertex[0]);
     } // if

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh	2009-05-26 21:59:31 UTC (rev 15059)
@@ -85,13 +85,10 @@
    * @param normalDir General preferred direction for fault normal
    *   (used to pick which of two possible normal directions for
    *   interface; only applies to fault surfaces in a 3-D domain).
-   * @param matDB Database of bulk elastic properties for fault region
-   *   (used to improve conditioning of Jacobian matrix)
    */
   void initialize(const topology::Mesh& mesh,
 		  const double upDir[3],
-		  const double normalDir[3],
-		  spatialdata::spatialdb::SpatialDB* matDB);
+		  const double normalDir[3]);
 
   /** Integrate contributions to residual term (r) for operator that
    * require assembly across processors.
@@ -188,15 +185,6 @@
   void _calcOrientation(const double upDir[3],
 			const double normalDir[3]);
 
-  /** Calculate conditioning field.
-   *
-   * @param cs Coordinate system for mesh
-   * @param matDB Database of bulk elastic properties for fault region
-   *   (used to improve conditioning of Jacobian matrix)
-   */
-  void _calcConditioning(const spatialdata::geocoords::CoordSys* cs,
-			 spatialdata::spatialdb::SpatialDB* matDB);
-
   /// Calculate fault area field.
   void _calcArea(void);
 

Modified: short/3D/PyLith/trunk/modulesrc/faults/Fault.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/Fault.i	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/modulesrc/faults/Fault.i	2009-05-26 21:59:31 UTC (rev 15059)
@@ -74,14 +74,11 @@
        * @param normalDir General preferred direction for fault normal
        *   (used to pick which of two possible normal directions for
        *   interface; only applies to fault surfaces in a 3-D domain).
-       * @param matDB Database of bulk elastic properties for fault region
-       *   (used to improve conditioning of Jacobian matrix)
        */
       virtual
       void initialize(const pylith::topology::Mesh& mesh,
 		      const double upDir[3],
-		      const double normalDir[3],
-		      spatialdata::spatialdb::SpatialDB* matDB) = 0;
+		      const double normalDir[3]) = 0;
       
       /** Get mesh associated with fault fields.
        *

Modified: short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i	2009-05-26 21:59:31 UTC (rev 15059)
@@ -42,13 +42,10 @@
        * @param normalDir General preferred direction for fault normal
        *   (used to pick which of two possible normal directions for
        *   interface; only applies to fault surfaces in a 3-D domain).
-       * @param matDB Database of bulk elastic properties for fault region
-       *   (used to improve conditioning of Jacobian matrix)
        */
       void initialize(const pylith::topology::Mesh& mesh,
 		      const double upDir[3],
-		      const double normalDir[3],
-		      spatialdata::spatialdb::SpatialDB* matDB);
+		      const double normalDir[3]);
 
       /** Integrate contribution of cohesive cells to residual term.
        *

Modified: short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveKin.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveKin.i	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveKin.i	2009-05-26 21:59:31 UTC (rev 15059)
@@ -58,13 +58,10 @@
        * @param normalDir General preferred direction for fault normal
        *   (used to pick which of two possible normal directions for
        *   interface; only applies to fault surfaces in a 3-D domain).
-       * @param matDB Database of bulk elastic properties for fault region
-       *   (used to improve conditioning of Jacobian matrix)
        */
       void initialize(const pylith::topology::Mesh& mesh,
 		      const double upDir[3],
-		      const double normalDir[3],
-		      spatialdata::spatialdb::SpatialDB* matDB);
+		      const double normalDir[3]);
       
       /** Integrate contributions to residual term (r) for operator that
        * require assembly across processors.

Modified: short/3D/PyLith/trunk/pylith/faults/Fault.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/Fault.py	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/pylith/faults/Fault.py	2009-05-26 21:59:31 UTC (rev 15059)
@@ -59,8 +59,6 @@
   @li \b normal_dir General preferred direction for fault normal
     (used to pick which of two possible normal directions for
     interface; only applies to fault surfaces in a 3-D domain).
-  @li \b mat_db Spatial database for bulk material properties
-    (used in improving conditioning of Jacobian matrix).
   
   \b Facilities
   @li \b quadrature Quadrature object for numerical integration
@@ -96,11 +94,6 @@
   faultQuadrature = pyre.inventory.facility("quadrature", factory=SubMeshQuadrature)
   faultQuadrature.meta['tip'] = "Quadrature object for numerical integration."
   
-  from spatialdata.spatialdb.SimpleDB import SimpleDB
-  matDB = pyre.inventory.facility("mat_db", family="spatial_database",
-                                  factory=SimpleDB)
-  matDB.meta['tip'] = "Spatial database for bulk material properties " \
-      "(used in improving conditioning of Jacobian matrix)."
 
   # PUBLIC METHODS /////////////////////////////////////////////////////
 
@@ -159,7 +152,7 @@
 
     self.faultQuadrature.initialize()
     ModuleFault.initialize(self, 
-                           self.mesh, self.upDir, self.normalDir, self.matDB)
+                           self.mesh, self.upDir, self.normalDir)
 
     if None != self.output:
       self.output.initialize(normalizer, self.faultQuadrature)
@@ -217,7 +210,6 @@
     self.faultQuadrature = self.inventory.faultQuadrature
     self.upDir = map(float, self.inventory.upDir)
     self.normalDir = map(float, self.inventory.normalDir)
-    self.matDB = self.inventory.matDB
     ModuleFault.id(self, self.inventory.matId)
     ModuleFault.label(self, self.inventory.faultLabel)
     return

Modified: short/3D/PyLith/trunk/tests/1d/line2/dislocation.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line2/dislocation.cfg	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/tests/1d/line2/dislocation.cfg	2009-05-26 21:59:31 UTC (rev 15059)
@@ -75,7 +75,6 @@
 [dislocation.timedependent.interfaces.fault]
 label = fault
 quadrature.cell.shape = point
-mat_db.iohandler.filename = matprops.spatialdb
 
 [dislocation.timedependent.interfaces.fault.eq_srcs.rupture.slip_function]
 slip.iohandler.filename = finalslip.spatialdb

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2009-05-26 21:59:31 UTC (rev 15059)
@@ -204,23 +204,6 @@
     CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->area[iVertex], areaVertex[0],
 				 tolerance);
   } // for
-
-  // Check pseudoStiffness
-  const ALE::Obj<RealSection>& stiffnessSection =
-    fault._fields->get("pseudostiffness").section();
-  CPPUNIT_ASSERT(!stiffnessSection.isNull());
-  iVertex = 0;
-  for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
-       v_iter != verticesEnd;
-       ++v_iter, ++iVertex) {
-    const int fiberDim = stiffnessSection->getFiberDimension(*v_iter);
-    CPPUNIT_ASSERT_EQUAL(1, fiberDim);
-    const double* stiffnessVertex = stiffnessSection->restrictPoint(*v_iter);
-    CPPUNIT_ASSERT(0 != stiffnessVertex);
-    const double tolerance = 1.0e-06;
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->pseudoStiffness, stiffnessVertex[0],
-				 tolerance);
-  } // for
 } // testInitialize
 
 // ----------------------------------------------------------------------
@@ -279,10 +262,9 @@
       
       const bool isConstraint = _isConstraintVertex(*v_iter);
       if (!isConstraint) {
-	const double pseudoStiffness = _data->pseudoStiffness;
 	for (int i=0; i < fiberDimE; ++i) {
 	  const int index = iVertex*spaceDim+i;
-	  const double valE = valsE[index] * pseudoStiffness;
+	  const double valE = valsE[index];
 	  if (fabs(valE) > tolerance)
 	    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
 	  else
@@ -318,10 +300,9 @@
       
       const bool isConstraint = _isConstraintVertex(*v_iter);
       if (!isConstraint) {
-	const double pseudoStiffness = _data->pseudoStiffness;
 	for (int i=0; i < fiberDimE; ++i) {
 	  const int index = iVertex*spaceDim+i;
-	  const double valE = valsE[index] * pseudoStiffness;
+	  const double valE = valsE[index];
 	  if (fabs(valE) > tolerance)
 	    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
 	  else
@@ -474,10 +455,9 @@
 	for (int i=0; i < fiberDimE; ++i)
 	  CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
       } else {
-	const double pseudoStiffness = _data->pseudoStiffness;
 	for (int i=0; i < fiberDimE; ++i) {
 	  const int index = iVertex*spaceDim+i;
-	  const double valE = valsE[index] * pseudoStiffness;
+	  const double valE = valsE[index];
 	  if (fabs(valE) > tolerance)
 	    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
 	  else
@@ -513,10 +493,9 @@
 	for (int i=0; i < fiberDimE; ++i)
 	  CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
       } else {
-	const double pseudoStiffness = _data->pseudoStiffness;
 	for (int i=0; i < fiberDimE; ++i) {
 	  const int index = iVertex*spaceDim+i;
-	  const double valE = valsE[index] * pseudoStiffness;
+	  const double valE = valsE[index];
 	  if (fabs(valE) > tolerance)
 	    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
 	  else
@@ -592,9 +571,7 @@
   for (int iRow=0; iRow < nrows; ++iRow)
     for (int iCol=0; iCol < ncols; ++iCol) {
       const int index = ncols*iRow+iCol;
-      const double pseudoStiffness = 
-	(iRow > iCol) ? 1.0 : _data->pseudoStiffness;
-      const double valE = valsE[index] * pseudoStiffness;
+      const double valE = valsE[index];
 #if 0 // DEBUGGING
       if (fabs(valE-vals[index]) > tolerance)
 	std::cout << "ERROR: iRow: " << iRow << ", iCol: " << iCol
@@ -797,7 +774,7 @@
     const double* dispVertex = dispSection->restrictPoint(meshVertex);
     CPPUNIT_ASSERT(0 != dispVertex);
 
-    const double scale = _data->pseudoStiffness / _data->area[iVertex];
+    const double scale = 1.0 / _data->area[iVertex];
     for (int iDim=0; iDim < spaceDim; ++iDim) {
       const double tractionE = dispVertex[iDim] * scale;
       if (tractionE > 1.0) 
@@ -886,13 +863,8 @@
   const double upDir[] = { 0.0, 0.0, 1.0 };
   const double normalDir[] = { 1.0, 0.0, 0.0 };
   
-  spatialdata::spatialdb::SimpleDB dbMatProp("material properties");
-  spatialdata::spatialdb::SimpleIOAscii ioMatProp;
-  ioMatProp.filename(_data->matPropsFilename);
-  dbMatProp.ioHandler(&ioMatProp);
+  fault->initialize(*mesh, upDir, normalDir); 
   
-  fault->initialize(*mesh, upDir, normalDir, &dbMatProp); 
-  
   delete[] sources; sources = 0;
   for (int i=0; i < nsrcs; ++i)
     delete[] names[i];

Modified: short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py	2009-05-26 21:06:59 UTC (rev 15058)
+++ short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py	2009-05-26 21:59:31 UTC (rev 15059)
@@ -310,14 +310,6 @@
     slipfn.inventory.dbSlipTime = dbSlipTime
     slipfn._configure()
 
-    ioMatDB = SimpleIOAscii()
-    ioMatDB.inventory.filename = "data/bulkprops_2d.spatialdb"
-    ioMatDB._configure()
-    dbMat = SimpleDB()
-    dbMat.inventory.iohandler = ioMatDB
-    dbMat.inventory.label = "bulk properties"
-    dbMat._configure()
-    
     # Setup fault
     fault = FaultCohesiveKin()
     fault.inventory.output.inventory.writer._configure()
@@ -327,7 +319,6 @@
     fault.inventory.upDir = [0, 0, 1]
     fault.inventory.normalDir = [1, 0, 0]
     fault.inventory.faultQuadrature = quadrature
-    fault.inventory.matDB = dbMat
     fault._configure()
     eqsrc = fault.eqsrcs.components()[0]
     eqsrc.inventory.originTime = 1.23*second



More information about the CIG-COMMITS mailing list