Thu Feb 23 17:11:00 PST 2012

Author: brad
Date: 2012-02-23 17:10:59 -0800 (Thu, 23 Feb 2012)
New Revision: 19668

Added:
   short/3D/PyLith/trunk/modulesrc/utils/constdefs.i
   short/3D/PyLith/trunk/unittests/pytests/utils/TestConstants.py
Removed:
   short/3D/PyLith/trunk/examples/3d/hex8/test.cfg
Modified:
   short/3D/PyLith/trunk/TODO
   short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
   short/3D/PyLith/trunk/examples/3d/hex8/step14.cfg
   short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
   short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
   short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/Integrator.cc
   short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc
   short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.hh
   short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc
   short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.hh
   short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc
   short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh
   short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.cc
   short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.hh
   short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.cc
   short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.hh
   short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.cc
   short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.hh
   short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticPlaneStrain.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticPlaneStress.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticStrain1D.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticStress1D.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellQpQsIsotropic3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/PowerLaw3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc
   short/3D/PyLith/trunk/libsrc/pylith/utils/constdefs.h
   short/3D/PyLith/trunk/modulesrc/friction/FrictionModel.i
   short/3D/PyLith/trunk/modulesrc/friction/RateStateAgeing.i
   short/3D/PyLith/trunk/modulesrc/friction/SlipWeakening.i
   short/3D/PyLith/trunk/modulesrc/friction/SlipWeakeningTime.i
   short/3D/PyLith/trunk/modulesrc/friction/StaticFriction.i
   short/3D/PyLith/trunk/modulesrc/friction/TimeWeakening.i
   short/3D/PyLith/trunk/modulesrc/utils/Makefile.am
   short/3D/PyLith/trunk/modulesrc/utils/utils.i
   short/3D/PyLith/trunk/tests/2d/frictionslide/plot_friction.py
   short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg
   short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
   short/3D/PyLith/trunk/unittests/libtests/friction/data/StaticFrictionData.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
   short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py
   short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc
   short/3D/PyLith/trunk/unittests/pytests/bc/TestAbsorbingDampers.py
   short/3D/PyLith/trunk/unittests/pytests/bc/TestNeumann.py
   short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveDyn.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/TestElasticityExplicitLgDeform.py
   short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py
   short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicitLgDeform.py
   short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticIsotropic3D.py
   short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticPlaneStrain.py
   short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticPlaneStress.py
   short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticStrain1D.py
   short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticStress1D.py
   short/3D/PyLith/trunk/unittests/pytests/topology/TestRefineUniform.py
   short/3D/PyLith/trunk/unittests/pytests/utils/Makefile.am
   short/3D/PyLith/trunk/unittests/pytests/utils/testutils.py
Log:
Merge from stable.

Modified: short/3D/PyLith/trunk/TODO
===================================================================

--- short/3D/PyLith/trunk/TODO	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/TODO	2012-02-24 01:10:59 UTC (rev 19668)
@@ -1,5 +1,5 @@
 ======================================================================
-CURRENT ISSUES/PRIORITIES (1.7.0)
+CURRENT ISSUES/PRIORITIES
 ======================================================================
 
 * --BUGS--
@@ -47,10 +47,6 @@
 
 * GenMaxwellQpQs [BRAD]
 
-  2-D and 3-D versions
-  Need to redo Maxwell time calculation. Use ratio.
-  Create benchmark for test (compare against fk)
-
 * Cleanup
 
     Add elasticPrestep() to Formulation (called from Problem)

Modified: short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg	2012-02-24 01:10:59 UTC (rev 19668)
@@ -74,7 +74,6 @@
 # PETSc
 # ----------------------------------------------------------------------
 # Set the solver options.
-
 [pylithapp.petsc]
 
 # Preconditioner settings.

Modified: short/3D/PyLith/trunk/examples/3d/hex8/step14.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/step14.cfg	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/examples/3d/hex8/step14.cfg	2012-02-24 01:10:59 UTC (rev 19668)
@@ -8,10 +8,13 @@
 #
 # This problem demonstrates the use of rate-and-state friction for a
 # quasi-static problem, using the aging law for evolution of the state
-# variable.
-# The problem is similar to example 13, except that a different friction
-# model is used. The model is run for 200 years. Slip begins to occur
-# at about 45 years, and continues in each step after that.
+# variable. The rate and state friction parameters result in
+# stick-slip behavior. 
+#
+# The problem is similar to example 13, except that a different
+# friction model is used. The model is run for 300 years with a time
+# step of 2.0 years. A smaller time step would better resolve the slip
+# time histories but increases the runtime.
 
 # ----------------------------------------------------------------------
 # RUNNING THE SIMULATON
@@ -49,11 +52,11 @@
 # solver.
 solver = pylith.problems.SolverNonlinear
 
-# Change the total simulation time to 200 years, and use a constant time
-# step size of 5 years.
+# Change the total simulation time to 300 years, and use a constant time
+# step size of 2 years.
 [pylithapp.timedependent.implicit.time_step]
-total_time = 200.0*year
-dt = 5.0*year
+total_time = 300.0*year
+dt = 2.0*year
 
 # ----------------------------------------------------------------------
 # boundary conditions
@@ -85,7 +88,7 @@
 db_initial = spatialdata.spatialdb.UniformDB
 db_initial.label = Dirichlet BC on +x
 db_initial.values = [displacement-x,displacement-y]
-db_initial.data = [-1.0*m,0.0*m]
+db_initial.data = [-2.0*m,0.0*m]
 
 db_rate = spatialdata.spatialdb.UniformDB
 db_rate.label = Dirichlet rate BC on +x
@@ -118,11 +121,12 @@
 [pylithapp.timedependent.interfaces.fault]
 # The label corresponds to the name of the nodeset in CUBIT.
 label = fault
-zero_tolerance = 1.0e-12
+zero_tolerance = 1.0e-11
 
 # Use the rate-and-state aging friction model.
 friction = pylith.friction.RateStateAgeing
 friction.label = Rate and state
+friction.min_slip_rate = 1.0e-9
 
 # We must define the quadrature information for fault cells.
 # The fault cells are 2D (surface).
@@ -131,22 +135,23 @@
 
 # Set rate-and-state parameters using a uniform DB. Set the parameters as
 # follows:
-# reference coefficient of friction: 0.6
-# reference slip rate: 1.0e-06 m/s
-# slip-weakening parameter: 0.037 m
-# a: 0.0125
-# b: 0.0172
+# reference coefficient of friction: 0.4
+# reference slip rate: 1.0e-11 m/s
+# characteristic slip distance: 0.05 m
+# a: 0.002
+# b: 0.08
 # cohesion: 0 Pa
 friction.db_properties = spatialdata.spatialdb.UniformDB
 friction.db_properties.label = Rate Stete Ageing
 friction.db_properties.values = [reference-friction-coefficient,reference-slip-rate,characteristic-slip-distance,constitutive-parameter-a,constitutive-parameter-b,cohesion]
-friction.db_properties.data = [0.6,1.0e-6*m/s,0.0370*m,0.0125,0.0172,0.0*Pa]
+friction.db_properties.data = [0.4,2.0e-11*m/s,0.05*m,0.002,0.08,0.0*Pa]
 
 # Set spatial database for the initial value of the state variable.
 friction.db_initial_state = spatialdata.spatialdb.UniformDB
 friction.db_initial_state.label = Rate State Ageing State
 friction.db_initial_state.values = [state-variable]
-friction.db_initial_state.data = [92.7*s]
+# theta_ss = characteristic_slip_dist / reference_slip_rate
+friction.db_initial_state.data = [2.5e+9*s]
 
 # ----------------------------------------------------------------------
 # PETSc settings
@@ -159,7 +164,7 @@
 # fault constitutive model.
 friction_pc_type = asm
 friction_sub_pc_factor_shift_type = nonzero
-friction_ksp_max_it = 25
+friction_ksp_max_it = 50
 friction_ksp_gmres_restart = 30
 # Uncomment to view details of friction sensitivity solve.
 #friction_ksp_monitor = true
@@ -167,13 +172,13 @@
 friction_ksp_converged_reason = true
 
 # Reduce convergence tolerances.
-ksp_rtol = 1.0e-13
-ksp_atol = 1.0e-15
+ksp_rtol = 1.0e-16
+ksp_atol = 1.0e-12
 
-snes_rtol = 1.0e-12
-snes_atol = 1.0e-14
+snes_rtol = 1.0e-14
+snes_atol = 1.0e-10
+snes_max_it = 200
 
-
 # ----------------------------------------------------------------------
 # output
 # ----------------------------------------------------------------------

Deleted: short/3D/PyLith/trunk/examples/3d/hex8/test.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/test.cfg	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/examples/3d/hex8/test.cfg	2012-02-24 01:10:59 UTC (rev 19668)
@@ -1,215 +0,0 @@
-# -*- Python -*-
-[pylithapp]
-
-# ----------------------------------------------------------------------
-# PROBLEM DESCRIPTION
-# ----------------------------------------------------------------------
-
-#
-# This problem demonstrates the use of rate-and-state friction for a
-# quasi-static problem, using the aging law for evolution of the state
-# variable.
-#
-# The problem is similar to example 13, except that a different
-# friction model is used. The axial deformation and time step have
-# also been changed. The model is run for 200 years.
-
-# ----------------------------------------------------------------------
-# RUNNING THE SIMULATON
-# ----------------------------------------------------------------------
-
-# This is not a self-contained simulation configuration file. This
-# file only specifies parameters specific to tutorial step14.
-# The general parameters are specificed in the pylithapp.cfg
-# file which PyLith reads by default.
-#
-# To run the simulation:
-# pylith step14.cfg
-#
-# Output will be directed to the directory output.
-
-# ----------------------------------------------------------------------
-# problem
-# ----------------------------------------------------------------------
-[pylithapp.timedependent]
-# Set bc to an array of 3 boundary conditions: 'x_pos', 'x_neg', and 'z_neg'.
-bc = [x_pos,x_neg,z_neg]
-
-# Set interfaces to an array of 1 fault: 'fault'.
-interfaces = [fault]
-
-[pylithapp.timedependent.implicit]
-# Set the output to an array of 2 output managers.
-# We will output the solution over the domain and the ground surface.
-output = [domain,subdomain]
-
-# Set subdomain component to OutputSolnSubset (boundary of the domain).
-output.subdomain = pylith.meshio.OutputSolnSubset
-
-# Fault friction is a nonlinear problem so we need to use the nonlinear
-# solver.
-solver = pylith.problems.SolverNonlinear
-
-# Change the total simulation time to 200 years, and use a constant time
-# step size of 5 years.
-[pylithapp.timedependent.implicit.time_step]
-total_time = 200.0*year
-dt = 0.5*year
-
-# ----------------------------------------------------------------------
-# boundary conditions
-# ----------------------------------------------------------------------
-# Set the parameters for boundary conditions applied on the
-# +x, -x, and -z faces of the box.
-#
-# On the -x and +x faces, we fix the x degrees of freedom and apply
-# velocities in the y-direction. We fix the z degree of freedom on the
-# bottom (-z) face.
-#
-# We use a UniformDB to apply the displacements and velocities, while
-# retaining the default ZeroDispDB for zero x-displacements on the -x
-# face and and zero z-displacements on -z.
-#
-# Note that since the fault cuts through the base of the model (z_neg),
-# we can only constrain the portion of the bottom boundary that does not
-# include the fault. A nodeset named 'face_zneg_nofault' has been defined
-# in Cubit for this purpose.
-#
-
-# The label corresponds to the name of the nodeset in CUBIT.
-
-# +x face -- Dirichlet
-[pylithapp.timedependent.bc.x_pos]
-bc_dof = [0,1]
-label = face_xpos
-
-db_initial = spatialdata.spatialdb.UniformDB
-db_initial.label = Dirichlet BC on +x
-db_initial.values = [displacement-x,displacement-y]
-db_initial.data = [-2.0*m,0.0*m]
-
-db_rate = spatialdata.spatialdb.UniformDB
-db_rate.label = Dirichlet rate BC on +x
-db_rate.values = [displacement-rate-x,displacement-rate-y,rate-start-time]
-db_rate.data = [0.0*cm/year,1.0*cm/year,0.0*year]
-
-# -x face
-[pylithapp.timedependent.bc.x_neg]
-bc_dof = [0, 1]
-label = face_xneg
-db_initial.label = Dirichlet BC on -x
-db_rate = spatialdata.spatialdb.UniformDB
-db_rate.label = Dirichlet rate BC on -x
-db_rate.values = [displacement-rate-x,displacement-rate-y,rate-start-time]
-db_rate.data = [0.0*cm/year,-1.0*cm/year,0.0*year]
-
-# -z face
-[pylithapp.timedependent.bc.z_neg]
-bc_dof = [2]
-label = face_zneg_nofault
-db_initial.label = Dirichlet BC on -z
-
-# ----------------------------------------------------------------------
-# faults
-# ----------------------------------------------------------------------
-[pylithapp.timedependent.interfaces]
-# Change fault to dynamic fault interface.
-fault = pylith.faults.FaultCohesiveDyn
-
-[pylithapp.timedependent.interfaces.fault]
-# The label corresponds to the name of the nodeset in CUBIT.
-label = fault
-zero_tolerance = 1.0e-12
-
-# Use the rate-and-state aging friction model.
-friction = pylith.friction.RateStateAgeing
-friction.label = Rate and state
-friction.min_slip_rate = 1.0e-9
-
-# We must define the quadrature information for fault cells.
-# The fault cells are 2D (surface).
-quadrature.cell = pylith.feassemble.FIATLagrange
-quadrature.cell.dimension = 2
-
-# Set rate-and-state parameters using a uniform DB. Set the parameters as
-# follows:
-# reference coefficient of friction: 0.4
-# reference slip rate: 1.0e-9 m/s
-# characteristic slip distance: 0.05 m
-# a: 0.01
-# b: 0.08
-# cohesion: 0 Pa
-friction.db_properties = spatialdata.spatialdb.UniformDB
-friction.db_properties.label = Rate Stete Ageing
-friction.db_properties.values = [reference-friction-coefficient,reference-slip-rate,characteristic-slip-distance,constitutive-parameter-a,constitutive-parameter-b,cohesion]
-friction.db_properties.data = [0.4, 1.0e-7*m/s, 2.0*m, 0.01, 0.02, 0.0*Pa]
-
-# Set spatial database for the initial value of the state variable.
-friction.db_initial_state = spatialdata.spatialdb.UniformDB
-friction.db_initial_state.label = Rate State Ageing State
-friction.db_initial_state.values = [state-variable]
-# theta_ss = characteristic_slip_dist / reference_slip_rate
-friction.db_initial_state.data = [5.0e+07*s]
-
-# ----------------------------------------------------------------------
-# PETSc settings
-# ----------------------------------------------------------------------
-# NOTE: There are additional settings specific to fault friction.
-[pylithapp.petsc]
-
-# Friction sensitivity solve used to compute the increment in slip
-# associated with changes in the Lagrange multiplier imposed by the
-# fault constitutive model.
-friction_pc_type = asm
-friction_sub_pc_factor_shift_type = nonzero
-friction_ksp_max_it = 50
-friction_ksp_gmres_restart = 30
-# Uncomment to view details of friction sensitivity solve.
-#friction_ksp_monitor = true
-#friction_ksp_view = true
-friction_ksp_converged_reason = true
-
-# Reduce convergence tolerances.
-ksp_rtol = 1.0e-16
-ksp_atol = 1.0e-15
-
-snes_rtol = 1.0e-14
-snes_atol = 1.0e-12
-
-snes_max_it = 200
-info = 1
-
-# ----------------------------------------------------------------------
-# output
-# ----------------------------------------------------------------------
-[pylithapp.problem.formulation.output.domain]
-output_freq = time_step
-time_step = 10.0*year
-writer = pylith.meshio.DataWriterHDF5Mesh
-writer.filename = output/test.h5
-
-[pylithapp.problem.formulation.output.subdomain]
-label = face_zpos
-writer = pylith.meshio.DataWriterHDF5SubMesh
-writer.filename = output/test-groundsurf.h5
-
-[pylithapp.problem.interfaces.fault.output]
-writer = pylith.meshio.DataWriterHDF5SubSubMesh
-writer.filename = output/test-fault.h5
-vertex_data_fields = [slip, slip_rate, traction, state_variable]
-
-[pylithapp.timedependent.materials.upper_crust.output]
-cell_filter = pylith.meshio.CellFilterAvgMesh
-output_freq = time_step
-time_step = 20.0*year
-writer = pylith.meshio.DataWriterHDF5Mesh
-writer.filename = output/test-upper_crust.h5
-
-# Give basename for VTK output of lower_crust state variables.
-[pylithapp.timedependent.materials.lower_crust.output]
-# Average values over quadrature points.
-cell_filter = pylith.meshio.CellFilterAvgMesh
-output_freq = time_step
-time_step = 20.0*year
-writer = pylith.meshio.DataWriterHDF5Mesh
-writer.filename = output/test-lower_crust.h5

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -338,6 +338,13 @@
       residualVertexN = areaVertex * (dispTpdtVertexL - initialTractionsVertex);
 
     } else { // opening, normal traction should be zero
+      if (fabs(tractionNormal) > _zeroTolerance) {
+	std::cerr << "WARNING! Fault opening with nonzero traction."
+		  << ", v_fault: " << v_fault
+		  << ", opening: " << slipNormal
+		  << ", normal traction: " << tractionNormal
+		  << std::endl;
+      } // if
       assert(fabs(tractionNormal) < _zeroTolerance);
     }  // if/else
     residualVertexP = -residualVertexN;
@@ -356,12 +363,6 @@
 	   residualSection->getFiberDimension(v_positive));
     residualSection->updateAddPoint(v_positive, &residualVertexP[0]);
 
-#if 0 // TEST
-    assert(residualVertexL.size() == 
-	   residualSection->getFiberDimension(v_lagrange));
-    residualSection->updateAddPoint(v_lagrange, &residualVertexL[0]);
-#endif
-
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(updateEvent);
 #endif
@@ -514,18 +515,11 @@
      const scalar_array&,
      const scalar_array&,
      const bool);
-  /// Member prototype for _constrainSolnSpaceImproveXD()
-  typedef void (pylith::faults::FaultCohesiveDyn::*constrainSolnSpaceImprove_fn_type)
-    (scalar_array*,
-     scalar_array*,
-     const scalar_array&,
-     const scalar_array&,
-     const scalar_array&);
 
-  assert(0 != fields);
-  assert(0 != _quadrature);
-  assert(0 != _fields);
-  assert(0 != _friction);
+  assert(fields);
+  assert(_quadrature);
+  assert(_fields);
+  assert(_friction);
 
   _sensitivitySetup(jacobian);
 
@@ -571,25 +565,18 @@
   assert(!dLagrangeTpdtSection.isNull());
 
   constrainSolnSpace_fn_type constrainSolnSpaceFn;
-  constrainSolnSpaceImprove_fn_type constrainSolnSpaceImproveFn;
   switch (spaceDim) { // switch
   case 1:
     constrainSolnSpaceFn = 
       &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D;
-    constrainSolnSpaceImproveFn = 
-      &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove1D;
     break;
   case 2: 
     constrainSolnSpaceFn = 
       &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D;
-    constrainSolnSpaceImproveFn = 
-      &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove2D;
     break;
   case 3:
     constrainSolnSpaceFn = 
       &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D;
-    constrainSolnSpaceImproveFn = 
-      &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove3D;
     break;
   default :
     assert(0);
@@ -749,18 +736,9 @@
     std::cout << ",  tractionTpdtVertex: ";
     for (int iDim=0; iDim < spaceDim; ++iDim)
       std::cout << "  " << tractionTpdtVertex[iDim];
-    std::cout << ",  lagrangeTVertex: ";
-    for (int iDim=0; iDim < spaceDim; ++iDim)
-      std::cout << "  " << lagrangeTVertex[iDim];
-    std::cout << ",  lagrangeTIncrVertex: ";
-    for (int iDim=0; iDim < spaceDim; ++iDim)
-      std::cout << "  " << lagrangeTIncrVertex[iDim];
     std::cout << ",  dTractionTpdtVertex: ";
     for (int iDim=0; iDim < spaceDim; ++iDim)
       std::cout << "  " << dTractionTpdtVertex[iDim];
-    std::cout << ",  dLagrangeTpdtVertex: ";
-    for (int iDim=0; iDim < spaceDim; ++iDim)
-      std::cout << "  " << dLagrangeTpdtVertex[iDim];
     std::cout << std::endl;
 #endif
      
@@ -798,30 +776,120 @@
   // change in Lagrange multipliers imposed by friction criterion.
 
   // Solve sensitivity problem for negative side of the fault.
-  bool negativeSide = true;
-  _sensitivityUpdateJacobian(negativeSide, jacobian, *fields);
-  _sensitivityReformResidual(negativeSide);
+  bool negativeSideFlag = true;
+  _sensitivityUpdateJacobian(negativeSideFlag, jacobian, *fields);
+  _sensitivityReformResidual(negativeSideFlag);
   _sensitivitySolve();
-  _sensitivityUpdateSoln(negativeSide);
+  _sensitivityUpdateSoln(negativeSideFlag);
 
   // Solve sensitivity problem for positive side of the fault.
-  negativeSide = false;
-  _sensitivityUpdateJacobian(negativeSide, jacobian, *fields);
-  _sensitivityReformResidual(negativeSide);
+  negativeSideFlag = false;
+  _sensitivityUpdateJacobian(negativeSideFlag, jacobian, *fields);
+  _sensitivityReformResidual(negativeSideFlag);
   _sensitivitySolve();
-  _sensitivityUpdateSoln(negativeSide);
+  _sensitivityUpdateSoln(negativeSideFlag);
 
   // Step 4: Update Lagrange multipliers and displacement fields based
   // on changes imposed by friction criterion in Step 2 (change in
   // Lagrange multipliers) and Step 3 (slip associated with change in
   // Lagrange multipliers).
+  //
+  // Use line search to find best update. This improves convergence
+  // because it accounts for feedback between the fault constitutive
+  // model and the deformation. We also search in log space because
+  // some fault constitutive models depend on the log of slip rate.
 
+  const PylithScalar residualTol = _zeroTolerance; // L2 misfit in tractions
+  const int maxIter = 16;
+  PylithScalar logAlphaL = log10(_zeroTolerance); // minimum step
+  PylithScalar logAlphaR = log10(1.0); // maximum step
+  PylithScalar logAlphaM = 0.5*(logAlphaL + logAlphaR);
+  PylithScalar logAlphaML = 0.5*(logAlphaL + logAlphaM);
+  PylithScalar logAlphaMR = 0.5*(logAlphaM + logAlphaR);
+  PylithScalar residualL = _constrainSolnSpaceNorm(pow(10.0, logAlphaL), t, fields);
+  PylithScalar residualML = _constrainSolnSpaceNorm(pow(10.0, logAlphaML), t, fields);
+  PylithScalar residualM = _constrainSolnSpaceNorm(pow(10.0, logAlphaM), t, fields);
+  PylithScalar residualMR = _constrainSolnSpaceNorm(pow(10.0, logAlphaMR), t, fields);
+  PylithScalar residualR = _constrainSolnSpaceNorm(pow(10.0, logAlphaR), t, fields);
+  for (int iter=0; iter < maxIter; ++iter) {
+    if (residualM < residualTol || residualR < residualTol)
+      // if residual is very small, we prefer the full step
+      break;
+
+#if 0
+    std::cout << "alphaL: " << pow(10.0, logAlphaL)
+	      << ", residuaL: " << residualL
+	      << ", alphaM: " << pow(10.0, logAlphaM)
+	      << ", residualM: " << residualM
+	      << ", alphaR: " << pow(10.0, logAlphaR)
+	      << ", residualR: " << residualR
+	      << std::endl;
+#endif
+
+    if (residualL < residualML && residualL < residualM && residualL < residualMR && residualL < residualR) {
+      logAlphaL = logAlphaL;
+      logAlphaR = logAlphaM;
+      residualL = residualL;
+      residualR = residualM;
+      residualM = residualML;
+    } else if (residualML <= residualL  && residualML < residualM && residualML < residualMR && residualML < residualR) {
+      logAlphaL = logAlphaL;
+      logAlphaR = logAlphaM;
+      residualL = residualL;
+      residualR = residualM;
+      residualM = residualML;
+    } else if (residualM <= residualL  && residualM <= residualML && residualM < residualMR && residualM < residualR) {
+      logAlphaL = logAlphaML;
+      logAlphaR = logAlphaMR;
+      residualL = residualML;
+      residualR = residualMR;
+      residualM = residualM;
+    } else if (residualMR <= residualL  && residualMR <= residualML && residualMR <= residualM && residualMR < residualR) {
+      logAlphaL = logAlphaM;
+      logAlphaR = logAlphaR;
+      residualL = residualM;
+      residualR = residualR;
+      residualM = residualMR;
+    } else if (residualR <= residualL  && residualR <= residualML && residualR <= residualM && residualR <= residualMR) {
+      logAlphaL = logAlphaM;
+      logAlphaR = logAlphaR;
+      residualL = residualM;
+      residualR = residualR;
+      residualM = residualMR;
+    } else {
+      assert(0);
+      throw std::logic_error("Unknown case in constrain solution space "
+			     "line search.");
+    } // if/else
+    logAlphaM = (logAlphaL + logAlphaR) / 2.0;
+    logAlphaML = (logAlphaL + logAlphaM) / 2.0;
+    logAlphaMR = (logAlphaM + logAlphaR) / 2.0;
+
+    residualML = _constrainSolnSpaceNorm(pow(10.0, logAlphaML), t, fields);
+    residualMR = _constrainSolnSpaceNorm(pow(10.0, logAlphaMR), t, fields);
+
+  } // for
+  // Account for possibility that end points have lowest residual
+  if (residualR <= residualM || residualR < residualTol) {
+    logAlphaM = logAlphaR;
+    residualM = residualR;
+  } else if (residualL < residualM) {
+    logAlphaM = logAlphaL;
+    residualM = residualL;
+  } // if/else
+  const PylithScalar alpha = pow(10.0, logAlphaM); // alphaM is our best guess
+#if 0 // DEBUGGING
+  std::cout << "ALPHA: " << alpha
+	    << ", residual: " << residualM
+	    << std::endl;
+#endif
+
   scalar_array slipTVertex(spaceDim);
   scalar_array dSlipTpdtVertex(spaceDim);
   scalar_array dispRelVertex(spaceDim);
+  const ALE::Obj<RealSection>& sensDispRelSection = _fields->get("sensitivity relative disp").section();
+  assert(!sensDispRelSection.isNull());
 
-  const ALE::Obj<RealSection>& sensDispRelSection =
-    _fields->get("sensitivity relative disp").section();
   for (int iVertex=0; iVertex < numVertices; ++iVertex) {
     const int v_fault = _cohesiveVertices[iVertex].fault;
     const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
@@ -832,14 +900,6 @@
     dLagrangeTpdtSection->restrictPoint(v_fault, &dLagrangeTpdtVertex[0],
 					dLagrangeTpdtVertex.size());
 
-    // If no change in the Lagrange multiplier computed from friction
-    // criterion, there are no updates, so continue.
-    PylithScalar dLagrangeTpdtMag = 0.0;
-    for (int i=0; i < spaceDim; ++i)
-      dLagrangeTpdtMag += dLagrangeTpdtVertex[i]*dLagrangeTpdtVertex[i];
-    if (0.0 == dLagrangeTpdtMag)
-      continue;
-
     // Get change in relative displacement from sensitivity solve.
     assert(spaceDim == sensDispRelSection->getFiberDimension(v_fault));
     const PylithScalar* sensDispRelVertex = 
@@ -887,11 +947,8 @@
       dispIncrSection->restrictPoint(v_lagrange);
     assert(lagrangeTIncrVertex);
 
-    // Step 4a: Prevent nonphysical trial solutions. The product of the
-    // normal traction and normal slip must be nonnegative (forbid
-    // interpenetration with tension or opening with compression).
-
-    // Compute slip, change in slip, and tractions in fault coordinates.
+    // Scale perturbation in relative displacements and change in
+    // Lagrange multipliers by alpha using only shear components.
     slipTVertex = 0.0;
     slipTpdtVertex = 0.0;
     dSlipTpdtVertex = 0.0;
@@ -905,25 +962,25 @@
 	  (dispTVertexP[jDim] - dispTVertexN[jDim] +
 	   dispIncrVertexP[jDim] - dispIncrVertexN[jDim]);
 	dSlipTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] * 
-	  sensDispRelVertex[jDim];
+	  alpha*sensDispRelVertex[jDim];
 	tractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
 	  (lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim]);
 	dTractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] * 
-	  dLagrangeTpdtVertex[jDim];
+	  alpha*dLagrangeTpdtVertex[jDim];
       } // for
     } // for
-    if (fabs(slipTpdtVertex[indexN]) < _zeroTolerance) {
-      slipTpdtVertex[indexN] = 0.0;
-    } // if
-    if (fabs(dSlipTpdtVertex[indexN]) < _zeroTolerance) {
-      dSlipTpdtVertex[indexN] = 0.0;
-    } // if
 
+    // FIRST, correct nonphysical trial solutions.
+    // Order of steps 5a-5c is important!
     if ((slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]) * 
 	(tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN])
 	< 0.0) {
+      // Step 5a: Prevent nonphysical trial solutions. The product of the
+      // normal traction and normal slip must be nonnegative (forbid
+      // interpenetration with tension or opening with compression).
+      
 #if 0 // DEBUGGING
-      std::cout << "STEP 4a CORRECTING NONPHYSICAL SLIP/TRACTIONS"
+      std::cout << "STEP 5a CORRECTING NONPHYSICAL SLIP/TRACTIONS"
 		<< ", v_fault: " << v_fault
 		<< ", slipNormal: " << slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]
 		<< ", tractionNormal: " << tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN]
@@ -941,33 +998,31 @@
 	dSlipTpdtVertex[indexN] = -slipTpdtVertex[indexN];
       } // if/else
 
-    } // if
-    // Do not allow fault interpenetration.
-    if (slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN] < 0.0) {
+    } else if (slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN] > _zeroTolerance) {
+      // Step 5b: Insure fault traction is zero when opening (if alpha=1
+      // this should be enforced already, but will not be properly
+      // enforced when alpha < 1).
+      for (int iDim=0; iDim < spaceDim; ++iDim) {
+	dTractionTpdtVertex[iDim] = -tractionTpdtVertex[iDim];
+      } // for
+      
+    } else if (slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN] < 0.0) {
+      // Step 5c: Prevent interpenetration.
 #if 0 // DEBUGGING
-      std::cout << "STEP 4a CORRECTING INTERPENETATION"
+      std::cout << "STEP 5b CORRECTING INTERPENETATION"
 		<< ", v_fault: " << v_fault
 		<< ", slipNormal: " << slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]
 		<< std::endl;
 #endif
       dSlipTpdtVertex[indexN] = -slipTpdtVertex[indexN];
-    } // if
+      
+    } // if/else      
 
-
-    // Improve estimate of slip and change in traction using dFriction/dD.
-    // Get friction properties and state variables.
-    _friction->retrievePropsStateVars(v_fault);
-
-    CALL_MEMBER_FN(*this,
-		   constrainSolnSpaceImproveFn)(&dTractionTpdtVertex, &dSlipTpdtVertex,
-						slipTVertex, slipTpdtVertex,
-						tractionTpdtVertex);
-
-#if 0 // OBSOLETE? Move to ImproveFn?
-    // Prevent over-correction in slip resulting in backslip
+#if 0 // UNNECESSARY????
+    // Step 5d: Prevent over-correction in shear slip resulting in backslip
     PylithScalar slipDot = 0.0;
     PylithScalar tractionSlipDot = 0.0;
-    for (int iDim=0; iDim < spaceDim-1; ++iDim)  { // :TODO: Optimize this
+    for (int iDim=0; iDim < indexN; ++iDim)  {
       // Compute dot product between slip and increment in slip (want positive)
       slipDot += 
 	(slipTpdtVertex[iDim] - slipTVertex[iDim]) * 
@@ -979,10 +1034,16 @@
     if (slipDot < 0.0 &&
 	sqrt(fabs(slipDot)) > _zeroTolerance && 
 	tractionSlipDot < 0.0) {
-      // Correct backslip
-      dTractionTpdtVertex *= 0.5; // Use bisection as guess for traction
-      for (int iDim=0; iDim < spaceDim-1; ++iDim) {
-	// Use bisection for slip
+#if 0 // DEBUGGING
+      std::cout << "STEP 5d CORRECTING BACKSLIP"
+		<< ", v_fault: " << v_fault
+		<< ", slipDot: " << slipDot
+		<< ", tractionSlipDot: " << tractionSlipDot
+		<< std::endl;
+#endif
+      // Correct backslip, use bisection as guess      
+      for (int iDim=0; iDim < indexN; ++iDim) {
+	dTractionTpdtVertex[iDim] *= 0.5;
 	dSlipTpdtVertex[iDim] = 0.5*(slipTVertex[iDim] - slipTpdtVertex[iDim]);
       } // for
     } // if/else
@@ -990,7 +1051,7 @@
     
     // Update current estimate of slip from t to t+dt.
     slipTpdtVertex += dSlipTpdtVertex;
-
+    
     // Compute relative displacement from slip.
     dispRelVertex = 0.0;
     dDispRelVertex = 0.0;
@@ -1044,7 +1105,6 @@
 
 #if 0 // DEBUGGING
   //dLagrangeTpdtSection->view("AFTER dLagrange");
-  dispRelSection->view("AFTER RELATIVE DISPLACEMENT");
   dispIncrSection->view("AFTER DISP INCR (t->t+dt)");
 #endif
 } // constrainSolnSpace
@@ -1393,8 +1453,6 @@
 #if 0 // DEBUGGING
   //dLagrangeTpdtSection->view("AFTER dLagrange");
   //dispIncrSection->view("AFTER DISP INCR (t->t+dt)");
-  dispRelSection->view("AFTER RELATIVE DISPLACEMENT");
-  //velRelSection->view("AFTER RELATIVE VELOCITY");
 #endif
 } // adjustSolnLumped
 
@@ -2182,8 +2240,13 @@
   scalar_array dispVertex(spaceDim);
   const ALE::Obj<RealSection>& solutionSection =
       _fields->get("sensitivity solution").section();
+  assert(!solutionSection.isNull());
   const ALE::Obj<RealSection>& dispRelSection =
     _fields->get("sensitivity relative disp").section();
+  assert(!dispRelSection.isNull());
+  const ALE::Obj<RealSection>& dLagrangeTpdtSection =
+      _fields->get("sensitivity dLagrange").section();
+  assert(!dLagrangeTpdtSection.isNull());
 
   const PylithScalar sign = (negativeSide) ? -1.0 : 1.0;
 
@@ -2193,6 +2256,20 @@
 
     solutionSection->restrictPoint(v_fault, &dispVertex[0], dispVertex.size());
 
+    // If no change in the Lagrange multiplier computed from friction
+    // criterion, there are no updates, so continue.
+    assert(spaceDim == dLagrangeTpdtSection->getFiberDimension(v_fault));
+    const PylithScalar* dLagrangeTpdtVertex = dLagrangeTpdtSection->restrictPoint(v_fault);
+    assert(dLagrangeTpdtVertex);
+    PylithScalar dLagrangeTpdtVertexMag = 0.0;
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
+      dLagrangeTpdtVertexMag += dLagrangeTpdtVertex[iDim]*dLagrangeTpdtVertex[iDim];
+    } // for
+    if (0.0 == dLagrangeTpdtVertexMag)
+      continue;
+
+    // Update relative displacements associated with sensitivity solve
+    // solution
     dispVertex *= sign;
 
     assert(dispVertex.size() == dispRelSection->getFiberDimension(v_fault));
@@ -2200,7 +2277,250 @@
   } // for
 } // _sensitivityUpdateSoln
 
+
 // ----------------------------------------------------------------------
+// Compute norm of residual associated with matching fault
+// constitutive model using update from sensitivity solve. We use
+// this in a line search to find a good update (required because
+// fault constitutive model may have a complex nonlinear feedback
+// with deformation).
+PylithScalar
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceNorm(const PylithScalar alpha,
+							  const PylithScalar t,
+							  topology::SolutionFields* const fields)
+{ // _constrainSolnSpaceNorm
+  /// Member prototype for _constrainSolnSpaceXD()
+  typedef void (pylith::faults::FaultCohesiveDyn::*constrainSolnSpace_fn_type)
+    (scalar_array*,
+     const PylithScalar,
+     const scalar_array&,
+     const scalar_array&,
+     const scalar_array&,
+     const bool);
+
+  // Update time step in friction (can vary).
+  _friction->timeStep(_dt);
+  const PylithScalar dt = _dt;
+
+  const int spaceDim = _quadrature->spaceDim();
+  const int indexN = spaceDim - 1;
+
+  constrainSolnSpace_fn_type constrainSolnSpaceFn;
+  switch (spaceDim) { // switch
+  case 1:
+    constrainSolnSpaceFn = 
+      &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D;
+    break;
+  case 2: 
+    constrainSolnSpaceFn = 
+      &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D;
+    break;
+  case 3:
+    constrainSolnSpaceFn = 
+      &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D;
+    break;
+  default :
+    assert(0);
+    throw std::logic_error("Unknown spatial dimension in "
+			   "FaultCohesiveDyn::constrainSolnSpace().");
+  } // switch
+
+  // Get sections
+  scalar_array slipTpdtVertex(spaceDim); // fault coordinates
+  scalar_array slipRateVertex(spaceDim); // fault coordinates
+  scalar_array tractionTpdtVertex(spaceDim); // fault coordinates
+  scalar_array tractionMisfitVertex(spaceDim); // fault coordinates
+
+  const ALE::Obj<RealSection>& orientationSection =
+      _fields->get("orientation").section();
+  assert(!orientationSection.isNull());
+
+  const ALE::Obj<RealSection>& dLagrangeTpdtSection =
+      _fields->get("sensitivity dLagrange").section();
+  assert(!dLagrangeTpdtSection.isNull());
+
+  const ALE::Obj<RealSection>& sensDispRelSection = 
+    _fields->get("sensitivity relative disp").section();
+  assert(!sensDispRelSection.isNull());
+
+  const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
+  assert(!dispTSection.isNull());
+
+  const ALE::Obj<RealSection>& dispIncrSection =
+      fields->get("dispIncr(t->t+dt)").section();
+  assert(!dispIncrSection.isNull());
+
+  bool isOpening = false;
+  PylithScalar norm2 = 0.0;
+  const int numVertices = _cohesiveVertices.size();
+  for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+    const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
+    const int v_fault = _cohesiveVertices[iVertex].fault;
+    const int v_negative = _cohesiveVertices[iVertex].negative;
+    const int v_positive = _cohesiveVertices[iVertex].positive;
+
+    // Get displacement values
+    assert(spaceDim == dispTSection->getFiberDimension(v_negative));
+    const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
+    assert(dispTVertexN);
+
+    assert(spaceDim == dispTSection->getFiberDimension(v_positive));
+    const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
+    assert(dispTVertexP);
+
+    // Get displacement increment values.
+    assert(spaceDim == dispIncrSection->getFiberDimension(v_negative));
+    const PylithScalar* dispIncrVertexN = 
+      dispIncrSection->restrictPoint(v_negative);
+    assert(dispIncrVertexN);
+
+    assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
+    const PylithScalar* dispIncrVertexP = 
+      dispIncrSection->restrictPoint(v_positive);
+    assert(dispIncrVertexP);
+
+    // Get orientation
+    assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+    const PylithScalar* orientationVertex = 
+      orientationSection->restrictPoint(v_fault);
+
+    // Get change in relative displacement from sensitivity solve.
+    assert(spaceDim == sensDispRelSection->getFiberDimension(v_fault));
+    const PylithScalar* dDispRelVertex = 
+      sensDispRelSection->restrictPoint(v_fault);
+    assert(dDispRelVertex);
+
+    // Get Lagrange multiplier values from disp(t), and dispIncr(t->t+dt)
+    assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
+    const PylithScalar* lagrangeTVertex = dispTSection->restrictPoint(v_lagrange);
+    assert(lagrangeTVertex);
+
+    assert(spaceDim == dispIncrSection->getFiberDimension(v_lagrange));
+    const PylithScalar* lagrangeTIncrVertex = 
+      dispIncrSection->restrictPoint(v_lagrange);
+    assert(lagrangeTIncrVertex);
+
+    assert(spaceDim == dLagrangeTpdtSection->getFiberDimension(v_fault));
+    const PylithScalar* dLagrangeTpdtVertex = 
+      dLagrangeTpdtSection->restrictPoint(v_fault);
+    assert(dLagrangeTpdtVertex);
+
+    // Compute slip, slip rate, and traction at time t+dt as part of
+    // line search.
+    slipTpdtVertex = 0.0;
+    slipRateVertex = 0.0;
+    tractionTpdtVertex = 0.0;
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
+      for (int jDim=0; jDim < spaceDim; ++jDim) {
+	slipTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+	  (dispTVertexP[jDim] + dispIncrVertexP[jDim]
+	   - dispTVertexN[jDim] - dispIncrVertexN[jDim] +
+	   alpha*dDispRelVertex[jDim]);
+	slipRateVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+	  (dispIncrVertexP[jDim] - dispIncrVertexN[jDim] +
+	   alpha*dDispRelVertex[jDim]) / dt;
+	tractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+	  (lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim] +
+	   alpha*dLagrangeTpdtVertex[jDim]);
+      } // for
+      if (fabs(slipRateVertex[iDim]) < _zeroTolerance) {
+	slipRateVertex[iDim] = 0.0;
+      } // if
+    } // for
+    if (fabs(slipTpdtVertex[indexN]) < _zeroTolerance) {
+      slipTpdtVertex[indexN] = 0.0;
+    } // if
+
+    // FIRST, correct nonphysical trial solutions.
+    // Order of steps a-c is important!
+
+    if (slipTpdtVertex[indexN]*tractionTpdtVertex[indexN] < 0.0) {
+      // Step a: Prevent nonphysical trial solutions. The product of the
+      // normal traction and normal slip must be nonnegative (forbid
+      // interpenetration with tension or opening with compression).
+      
+      // Don't know what behavior is appropriate so set smaller of
+      // traction and slip to zero (should be appropriate if problem
+      // is nondimensionalized correctly).
+      if (fabs(slipTpdtVertex[indexN]) > fabs(tractionTpdtVertex[indexN])) {
+	// fault opening is bigger, so force normal traction back to zero
+	tractionTpdtVertex[indexN] = 0.0;
+      } else {
+	// traction is bigger, so force fault opening back to zero
+	slipTpdtVertex[indexN] = 0.0;
+      } // if/else
+
+    } else if (slipTpdtVertex[indexN] > _zeroTolerance) {
+      // Step b: Insure fault traction is zero when opening (if
+      // alpha=1 this should be enforced already, but will not be
+      // properly enforced when alpha < 1).
+      
+      for (int iDim=0; iDim < spaceDim; ++iDim) {
+	tractionTpdtVertex[iDim] = 0.0;
+      } // for
+    } else if (slipTpdtVertex[indexN] < 0.0) {
+      // Step c: Prevent interpenetration.
+
+      slipTpdtVertex[indexN] = 0.0;
+    } // if
+
+    if (slipTpdtVertex[indexN] > _zeroTolerance) {
+      isOpening = true;
+    } // if
+
+    // Apply friction criterion to trial solution to get change in
+    // Lagrange multiplier (dLagrangeTpdtVertex) in fault coordinate
+    // system.
+    
+    // Get friction properties and state variables.
+    _friction->retrievePropsStateVars(v_fault);
+    
+    // Use fault constitutive model to compute traction associated with
+    // friction.
+    tractionMisfitVertex = 0.0;
+    const bool iterating = true; // Iterating to get friction
+    CALL_MEMBER_FN(*this,
+		   constrainSolnSpaceFn)(&tractionMisfitVertex, t,
+					 slipTpdtVertex, slipRateVertex,
+					 tractionTpdtVertex,
+					 iterating);
+
+#if 0 // DEBUGGING
+    std::cout << "alpha: " << alpha
+	      << ", v_fault: " << v_fault;
+    std::cout << ", misfit:";
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
+      std::cout << " " << tractionMisfitVertex[iDim];
+    } // for
+    std::cout << ", slip:";
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
+      std::cout << " " << slipTpdtVertex[iDim];
+    } // for
+    std::cout << ", traction:";
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
+      std::cout << " " << tractionTpdtVertex[iDim];
+    } // for
+    std::cout << ", dDispRel:";
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
+      std::cout << " " << dDispRelVertex[iDim];
+    } // for
+    std::cout << std::endl;
+#endif
+
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
+      norm2 += tractionMisfitVertex[iDim]*tractionMisfitVertex[iDim];
+    } // for
+  } // for
+
+  if (isOpening && alpha < 1.0) {
+    norm2 = PYLITH_MAXFLOAT;
+  } // if
+
+  return sqrt(norm2) / numVertices;
+} // _constrainSolnSpaceNorm
+
+
+// ----------------------------------------------------------------------
 // Constrain solution space in 1-D.
 void
 pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D(scalar_array* dTractionTpdt,
@@ -2338,320 +2658,4 @@
 } // _constrainSolnSpace3D
 
 
-// ----------------------------------------------------------------------
-// Constrain solution space in 1-D.
-void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove1D(
-	 scalar_array* dTractionTpdt,
-	 scalar_array* dSlipTpdt,
-         const scalar_array& slipT,
-         const scalar_array& slipTpdt,
-	 const scalar_array& tractionTpdt)
-{ // _constrainSolnSpaceImprove3D
-  assert(dTractionTpdt);
-  assert(dSlipTpdt);
-
-  // Improving slip estimate only applies in shear. Do nothing.
-
-  PetscLogFlops(0); // :TODO: Update this.
-} // _constrainSolnSpaceImprove1D
-
-
-// ----------------------------------------------------------------------
-// Constrain solution space in 3-D.
-void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove2D(
-	 scalar_array* dTractionTpdt,
-	 scalar_array* dSlipTpdt,
-         const scalar_array& slipT,
-         const scalar_array& slipTpdt,
-	 const scalar_array& tractionTpdt)
-{ // _constrainSolnSpaceImprove2D
-  assert(dTractionTpdt);
-  assert(dSlipTpdt);
-
-#if 0 // DEBUGGING
-  std::cout << "BEFORE improvement"
-	    << ", dTractionTpdt:";
-  for (int i=0; i < 2; ++i)
-    std::cout << "  " << (*dTractionTpdt)[i];
-  std::cout << ", dSlipTpdt:";
-  for (int i=0; i < 2; ++i)
-    std::cout << "  " << (*dSlipTpdt)[i];
-  std::cout << std::endl;
-#endif
-
-  // Compute magnitude of slip and slip rate (with current increment).
-  const PylithScalar slipShearMag = fabs(slipTpdt[0]);
-  const PylithScalar slipShearNewMag = fabs(slipTpdt[0]+(*dSlipTpdt)[0]);
-  const PylithScalar slipRateMag = fabs(slipTpdt[0]-slipT[0]) / _dt;
-  const PylithScalar dSlipShearNewMag = fabs((*dSlipTpdt)[0]);
-
-  const PylithScalar tractionNormal = tractionTpdt[1];
-
-  // Friction stress for old estimate of slip is tractionTpdt +
-  // dTractionTpdt.
-  const PylithScalar frictionStress = 
-    fabs(tractionTpdt[0]+(*dTractionTpdt)[0]);
-  const PylithScalar tractionShearMag = fabs(tractionTpdt[0]);
-  
-  if (fabs(slipTpdt[1] + (*dSlipTpdt)[1]) < _zeroTolerance &&
-      tractionNormal < -_zeroTolerance &&
-      dSlipShearNewMag > 0.0) {
-    // if in compression and no opening, and changing slip
-
-    // Calculate slope (Jacobian) of friction at slip before adding
-    // new increment.
-    const PylithScalar slopeF = 
-      _friction->calcFrictionSlope(slipShearMag, slipRateMag, tractionNormal);
-
-#if 0 // linear space
-    const PylithScalar slopeT = 
-      (frictionStress - tractionShearMag) / dSlipShearNewMag;
-
-    // Set adjustments to increments to original increments as default.
-    PylithScalar dSlipShearNew2Mag = dSlipShearNewMag;
-    PylithScalar dTractionShearNew2Mag = frictionStress - tractionShearMag;
-    if (slopeF > 0.0 && tractionShearMag > frictionStress) {
-      // Strengthening, so reduce increment in slip
-      dSlipShearNew2Mag = 
-	(tractionShearMag - frictionStress) / (slopeF-slopeT));
-      dTractionShearNew2Mag += slopeF * dSlipShearNew2Mag;
-    } // if
-    // Ignore other cases, because slip will exceed estimate based on
-    // elasticity
-
-#else // log space
-
-    const PylithScalar slopeT = (frictionStress - tractionShearMag) / 
-      (log(slipShearNewMag) - log(slipShearMag));
-    
-    // Set adjustments to increments to original increments as default.
-    PylithScalar dSlipShearNew2Mag = dSlipShearNewMag;
-    PylithScalar dTractionShearNew2Mag = frictionStress - tractionShearMag;
-    if (slopeF > 0.0 && tractionShearMag > frictionStress) {
-      // Strengthening, so reduce increment in slip
-      const PylithScalar slipShearMagEff = std::max(slipShearMag, _zeroTolerance);
-      dSlipShearNew2Mag = slipShearMagEff * 
-	(-1.0 + exp((tractionShearMag - frictionStress) / (slopeF-slopeT)));
-      dTractionShearNew2Mag += 
-	slopeF * log((slipShearMagEff + dSlipShearNew2Mag)/slipShearMagEff);
-
-    } // if
-    // Ignore other cases, because slip will exceed estimate based on
-    // elasticity
-#endif
-    
-    // Project slip and traction into vector components
-    (*dSlipTpdt)[0] *= dSlipShearNew2Mag / dSlipShearNewMag;
-
-    const PylithScalar dTractionTpdtMag = fabs((*dTractionTpdt)[0]);
-    assert(dTractionTpdtMag > 0.0);
-    (*dTractionTpdt)[0] *= fabs(dTractionShearNew2Mag) / dTractionTpdtMag;
-
-    // Prevent over-correction in slip resulting in backslip.
-    // Expect slip direction to match tractions
-
-    // Compute dot product between slip and increment in slip (want positive)
-    const PylithScalar slipDot = 
-      (slipTpdt[0] - slipT[0]) * (slipTpdt[0] + (*dSlipTpdt)[0] - slipT[0]);
-    // Compute dot product of traction and slip
-    const PylithScalar tractionSlipDot = 
-      (tractionTpdt[0] + (*dTractionTpdt)[0]) * (slipTpdt[0] + (*dSlipTpdt)[0]);
-    if (slipDot < 0.0 &&
-	sqrt(fabs(slipDot)) > _zeroTolerance && 
-	tractionSlipDot < 0.0) {
-      // Correct backslip
-      (*dTractionTpdt)[0] *= 0.5; // Use bisection as guess for traction
-      // Use bisection for slip
-#if 0 // linear space
-      (*dSlipTpdt)[0] = 0.5*(slipT[0] - slipTpdt[0]);
-#else // log space
-      assert(slipT[0] * slipTpdt[0] > 0.0);
-      (*dSlipTpdt)[0] *= sqrt(slipT[0] * slipTpdt[0]) / dSlipShearNew2Mag;
-#endif
-    } // if/else
-
-#if 0 // DEBUGGING
-  std::cout << "AFTER improvement"
-	    << ", dTractionTpdt:";
-  for (int i=0; i < 2; ++i)
-    std::cout << "  " << (*dTractionTpdt)[i];
-  std::cout << ", dSlipTpdt:";
-  for (int i=0; i < 2; ++i)
-    std::cout << "  " << (*dSlipTpdt)[i];
-  std::cout << ", frictionStress: " << frictionStress
-	    << ", tractionShearMag: " << tractionShearMag
-	    << ", slopeF: " << slopeF
-	    << ", slopeT: " << slopeT
-	    << std::endl;
-#endif
-
-  } // if
-
-  PetscLogFlops(0); // :TODO: Update this.
-
-
-} // _constrainSolnSpaceImprove2D
-
-
-// ----------------------------------------------------------------------
-// Constrain solution space in 3-D.
-void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove3D(
-	 scalar_array* dTractionTpdt,
-	 scalar_array* dSlipTpdt,
-         const scalar_array& slipT,
-         const scalar_array& slipTpdt,
-	 const scalar_array& tractionTpdt)
-{ // _constrainSolnSpaceImprove3D
-  assert(dTractionTpdt);
-  assert(dSlipTpdt);
- 
-#if 0 // DEBUGGING
-  std::cout << "BEFORE improvement"
-	    << ", dTractionTpdt:";
-  for (int i=0; i < 2; ++i)
-    std::cout << "  " << (*dTractionTpdt)[i];
-  std::cout << ", dSlipTpdt:";
-  for (int i=0; i < 2; ++i)
-    std::cout << "  " << (*dSlipTpdt)[i];
-  std::cout << std::endl;
-#endif
-
- // Compute magnitude of slip and slip rate (with current increment).
-  const PylithScalar slipShearMag = 
-    sqrt(slipTpdt[0]*slipTpdt[0] +
-	 slipTpdt[1]*slipTpdt[1]);
-  const PylithScalar slipShearNewMag = 
-    sqrt(pow(slipTpdt[0]+(*dSlipTpdt)[0], 2) +
-	 pow(slipTpdt[1]+(*dSlipTpdt)[1], 2));
-  const PylithScalar slipRateMag = 
-    sqrt(pow(slipTpdt[0]-slipT[0], 2) +
-	 pow(slipTpdt[1]-slipT[1], 2)) / _dt;
-  const PylithScalar dSlipShearNewMag = 
-    sqrt(pow((*dSlipTpdt)[0], 2) +
-	 pow((*dSlipTpdt)[1], 2));
-
-    // Friction stress for old estimate of slip is tractionTpdt +
-    // dTractionTpdt.
-    const PylithScalar frictionStress = 
-      sqrt(pow(tractionTpdt[0]+(*dTractionTpdt)[0], 2) +
-	   pow(tractionTpdt[1]+(*dTractionTpdt)[1], 2));
-    const PylithScalar tractionShearMag = 
-      sqrt(tractionTpdt[0]*tractionTpdt[0] +
-	   tractionTpdt[1]*tractionTpdt[1]);
-
-  const PylithScalar tractionNormal = tractionTpdt[2] + (*dTractionTpdt)[2];
-
-  if (fabs(slipTpdt[2] + (*dSlipTpdt)[2]) < _zeroTolerance &&
-      tractionNormal < -_zeroTolerance &&
-      dSlipShearNewMag > 0.0) {
-    // if in compression and no opening, and changing slip
-
-    // Calculate slope (Jacobian) of friction at slip before adding
-    // new increment.
-    const PylithScalar slopeF = 
-      _friction->calcFrictionSlope(slipShearMag, slipRateMag, tractionNormal);
-
-#if 0 // linear space
-    const PylithScalar slopeT = 
-      (frictionStress - tractionShearMag) / dSlipShearNewMag;
-
-    // Set adjustments to increments to original increments as default.
-    PylithScalar dSlipShearNew2Mag = dSlipShearNewMag;
-    PylithScalar dTractionShearNew2Mag = frictionStress - tractionShearMag;
-    if (slopeF > 0.0 && tractionShearMag > frictionStress) {
-      // Strengthening, so reduce increment in slip
-      dSlipShearNew2Mag = 
-	(tractionShearMag - frictionStress) / (slopeF-slopeT));
-      dTractionShearNew2Mag += slopeF * dSlipShearNew2Mag;
-    } // if
-    // Ignore other cases, because slip will exceed estimate based on
-    // elasticity
-
-#else // log space
-
-    const PylithScalar slopeT = (frictionStress - tractionShearMag) / 
-      (log(slipShearNewMag) - log(slipShearMag));
-    
-    // Set adjustments to increments to original increments as default.
-    PylithScalar dSlipShearNew2Mag = dSlipShearNewMag;
-    PylithScalar dTractionShearNew2Mag = frictionStress - tractionShearMag;
-    if (slopeF > 0.0 && tractionShearMag > frictionStress) {
-      // Strengthening, so reduce increment in slip
-      const PylithScalar slipShearMagEff = std::max(slipShearMag, _zeroTolerance);
-      dSlipShearNew2Mag = slipShearMagEff * 
-	(-1.0 + exp((tractionShearMag - frictionStress) / (slopeF-slopeT)));
-      dTractionShearNew2Mag += 
-	slopeF * log((slipShearMagEff + dSlipShearNew2Mag)/slipShearMagEff);
-
-    } // if
-    // Ignore other cases, because slip will exceed estimate based on
-    // elasticity
-#endif
-    
-    // Project slip and traction into vector components
-    // keeping same direction as original
-    (*dSlipTpdt)[0] *= dSlipShearNew2Mag / dSlipShearNewMag;
-    (*dSlipTpdt)[1] *= dSlipShearNew2Mag / dSlipShearNewMag;
-
-    const PylithScalar dTractionTpdtMag = 
-      sqrt(pow((*dTractionTpdt)[0], 2) +
-	   pow((*dTractionTpdt)[1], 2));
-    assert(dTractionTpdtMag > 0.0);
-    (*dTractionTpdt)[0] *= fabs(dTractionShearNew2Mag) / dTractionTpdtMag;
-    (*dTractionTpdt)[1] *= fabs(dTractionShearNew2Mag) / dTractionTpdtMag;
-
-    // Prevent over-correction in slip resulting in backslip.
-    // Expect slip direction to match tractions
-
-    // Compute dot product between slip and increment in slip (want positive)
-    const PylithScalar slipDot = 
-      (slipTpdt[0] - slipT[0]) * (slipTpdt[0] + (*dSlipTpdt)[0] - slipT[0]) +
-      (slipTpdt[1] - slipT[1]) * (slipTpdt[1] + (*dSlipTpdt)[1] - slipT[1]);
-    // Compute dot product of traction and slip
-    const PylithScalar tractionSlipDot = 
-      (tractionTpdt[0] + (*dTractionTpdt)[0]) * (slipTpdt[0] + (*dSlipTpdt)[0])+
-      (tractionTpdt[1] + (*dTractionTpdt)[1]) * (slipTpdt[1] + (*dSlipTpdt)[1]);
-    if (slipDot < 0.0 &&
-	sqrt(fabs(slipDot)) > _zeroTolerance && 
-	tractionSlipDot < 0.0) {
-      // Correct backslip
-      (*dTractionTpdt)[0] *= 0.5; // Use bisection as guess for traction
-      (*dTractionTpdt)[1] *= 0.5;
-      // Use bisection for slip
-#if 0 // linear space
-      (*dSlipTpdt)[0] = 0.5*(slipT[0] - slipTpdt[0]);
-      (*dSlipTpdt)[1] = 0.5*(slipT[1] - slipTpdt[1]);
-#else // log space
-      assert(slipT[0] * slipTpdt[0] > 0.0);
-      assert(slipT[1] * slipTpdt[1] > 0.0);
-
-      (*dSlipTpdt)[0] *= sqrt(slipT[0] * slipTpdt[0]) / dSlipShearNew2Mag; 
-      (*dSlipTpdt)[1] *= sqrt(slipT[1] * slipTpdt[1]) / dSlipShearNew2Mag;
-#endif
-    } // if
-
-#if 0 // DEBUGGING
-  std::cout << "AFTER improvement"
-	    << ", dTractionTpdt:";
-  for (int i=0; i < 2; ++i)
-    std::cout << "  " << (*dTractionTpdt)[i];
-  std::cout << ", dSlipTpdt:";
-  for (int i=0; i < 2; ++i)
-    std::cout << "  " << (*dSlipTpdt)[i];
-  std::cout << ", frictionStress: " << frictionStress
-	    << ", tractionShearMag: " << tractionShearMag
-	    << ", slopeF: " << slopeF
-	    << ", slopeT: " << slopeT
-	    << std::endl;
-#endif
-
-} // if
-
-  PetscLogFlops(0); // :TODO: Update this.
-} // _constrainSolnSpaceImprove3D
-
-
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh	2012-02-24 01:10:59 UTC (rev 19668)
@@ -201,6 +201,22 @@
    */
   void _sensitivityUpdateSoln(const bool negativeSide);
 
+  /** Compute norm of residual associated with matching fault
+   *  constitutive model using update from sensitivity solve. We use
+   *  this in a line search to find a good update (required because
+   *  fault constitutive model may have a complex nonlinear feedback
+   *  with deformation).
+   *
+   * @param alpha Current step in line search.
+   * @param t Current time.
+   * @param fields Solution fields.
+   *
+   * @returns L2 norm of residual.
+   */
+  PylithScalar _constrainSolnSpaceNorm(const PylithScalar alpha,
+				       const PylithScalar t,
+				       topology::SolutionFields* const fields);
+
   /** Constrain solution space in 1-D.
    *
    * @param dLagrangeTpdt Adjustment to Lagrange multiplier.
@@ -246,48 +262,6 @@
 			     const scalar_array& tractionTpdt,
 			     const bool iterating =true);
 
-  /** Improve slip estimate when constraining solution space in 1-D.
-   *
-   * @param dTractionTpdt Adjustment to fault traction vector.
-   * @param dSlipTpdt Adjustment to fault slip vector.
-   * @param slipT Fault slip vector at time t.
-   * @param slipTpdt Fault slip vector at time t+dt (without adjustment).
-   * @param tractionTpdt Fault traction vector (without adjustment).
-   */
-  void _constrainSolnSpaceImprove1D(double_array* dTractionTpdt,
-				    double_array* dSlipTpdt,
-				    const double_array& slipT,
-				    const double_array& slipTpdt,
-				    const double_array& tractionTpdt);
-
-  /** Improve slip estimate when constraining solution space in 2-D.
-   *
-   * @param dTractionTpdt Adjustment to fault traction vector.
-   * @param dSlipTpdt Adjustment to fault slip vector.
-   * @param slipT Fault slip vector at time t.
-   * @param slipTpdt Fault slip vector at time t+dt (without adjustment).
-   * @param tractionTpdt Fault traction vector (without adjustment).
-   */
-  void _constrainSolnSpaceImprove2D(double_array* dTractionTpdt,
-				    double_array* dSlipTpdt,
-				    const double_array& slipT,
-				    const double_array& slipTpdt,
-				    const double_array& tractionTpdt);
-
-  /** Improve slip estimate when constraining solution space in 3-D.
-   *
-   * @param dTractionTpdt Adjustment to fault traction vector.
-   * @param dSlipTpdt Adjustment to fault slip vector.
-   * @param slipT Fault slip vector at time t.
-   * @param slipTpdt Fault slip vector at time t+dt (without adjustment).
-   * @param tractionTpdt Fault traction vector (without adjustment).
-   */
-  void _constrainSolnSpaceImprove3D(double_array* dTractionTpdt,
-				    double_array* dSlipTpdt,
-				    const double_array& slipT,
-				    const double_array& slipTpdt,
-				    const double_array& tractionTpdt);
-
   // PRIVATE MEMBERS ////////////////////////////////////////////////////
 private :
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -1876,7 +1876,7 @@
     const ALE::Obj<RealSection>& partitionSection = partition.section();
     assert(!partitionSection.isNull());
     
-    const double rank = (double) partitionSection->commRank();
+    const PylithScalar rank = (double) partitionSection->commRank();
     // Loop over cells in fault mesh, set partition
     for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin; 
 	 c_iter != cellsEnd;

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Integrator.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Integrator.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -20,7 +20,7 @@
 
 #include "Quadrature.hh" // USES Quadrature
 #include "pylith/utils/EventLogger.hh" // USES EventLogger
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 #include "spatialdata/spatialdb/GravityField.hh" // USES GravityField
@@ -113,7 +113,7 @@
 pylith::feassemble::Integrator<quadrature_type>::stableTimeStep(const topology::Mesh& mesh)
 { // stableTimeStep
   // Assume any time step will work.
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 } // stableTimeStep
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -325,28 +325,6 @@
 } // calcFriction
 
 // ----------------------------------------------------------------------
-// Compute change in friction with change in slip.
-PylithScalar
-pylith::friction::FrictionModel::calcFrictionSlope(const PylithScalar slip,
-						   const PylithScalar slipRate,
-						   const PylithScalar normalTraction)
-{ // calcFrictionSlope
-  assert(_fieldsPropsStateVars);
-  
-  assert(_propsFiberDim+_varsFiberDim == _propsStateVarsVertex.size());
-  const PylithScalar* propertiesVertex = &_propsStateVarsVertex[0];
-  const PylithScalar* stateVarsVertex = (_varsFiberDim > 0) ?
-    &_propsStateVarsVertex[_propsFiberDim] : 0;
-  
-  const PylithScalar slope =
-    _calcFrictionSlope(slip, slipRate, normalTraction,
-		       propertiesVertex, _propsFiberDim,
-		       stateVarsVertex, _varsFiberDim);
-  
-  return slope;
-} // calcFrictionSlope
-
-// ----------------------------------------------------------------------
 // Update state variables (for next time step).
 void
 pylith::friction::FrictionModel::updateStateVars(const PylithScalar t,

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.hh	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.hh	2012-02-24 01:10:59 UTC (rev 19668)
@@ -168,21 +168,6 @@
 			    const PylithScalar slipRate,
 			    const PylithScalar normalTraction);
   
-  /** Compute change in friction for a change in slip (Jacobian).
-   *
-   * @pre Must call retrievePropsAndVars for cell before calling
-   * calcFriction().
-   *
-   * @param slip Current slip at location.
-   * @param slipRate Current slip rate at location.
-   * @param normalTraction Normal traction at location.
-   *
-   * @returns Change in friction for a chance in slip (dT/dD).
-   */
-  PylithScalar calcFrictionSlope(const PylithScalar slip,
-			   const PylithScalar slipRate,
-			   const PylithScalar normalTraction);
-  
   /** Compute friction at vertex.
    *
    * @pre Must call retrievePropsAndVars for cell before calling
@@ -282,27 +267,6 @@
 			     const PylithScalar* stateVars,
 			     const int numStateVars) = 0;
   
-  /** Compute change in friction for a change in slip (Jacobian).
-   *
-   * @param slip Current slip at location.
-   * @param slipRate Current slip rate at location.
-   * @param normalTraction Normal traction at location.
-   * @param properties Properties at location.
-   * @param numProperties Number of properties.
-   * @param stateVars State variables at location.
-   * @param numStateVars Number of state variables.
-   *
-   * @returns Change in friction for a chance in slip (dT/dD).
-   */
-  virtual
-  PylithScalar _calcFrictionSlope(const PylithScalar slip,
-			    const PylithScalar slipRate,
-			    const PylithScalar normalTraction,
-			    const PylithScalar* properties,
-			    const int numProperties,
-			    const PylithScalar* stateVars,
-			    const int numStateVars) = 0;
-  
   /** Update state variables (for next time step).
    *
    * @param t Time in simulation.

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "pylith/materials/Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -310,31 +310,6 @@
   if (normalTraction <= 0.0) {
     // if fault is in compression
 
-#if 0
-    // regularized rate and state equation
-    const PylithScalar f0 = properties[p_coef];
-
-    // Since regulatized friction -> 0 as slipRate -> 0, limit slip
-    // rate to some minimum value
-    const PylithScalar slipRateEff = std::max(_minSlipRate, slipRate);
-
-    const PylithScalar slipRate0 = properties[p_slipRate0];
-    const PylithScalar a = properties[p_a];
-
-    const PylithScalar theta = stateVars[s_state];
-
-    const PylithScalar L = properties[p_L];
-    const PylithScalar b = properties[p_b];
-    const PylithScalar bLnTerm = b * log(slipRate0 * theta / L);
-
-    const PylithScalar expTerm = exp((f0 + bLnTerm)/a);
-    const PylithScalar sinhArg = 0.5 * slipRateEff / slipRate0 * expTerm;
-
-    mu_f = a * asinh(sinhArg);
-    friction = -mu_f * normalTraction + properties[p_cohesion];
-
-#else
-
     const PylithScalar slipRateLinear = _minSlipRate;
 
     const PylithScalar f0 = properties[p_coef];
@@ -350,61 +325,18 @@
       mu_f = f0 + a*log(slipRateLinear / slipRate0) + b*log(slipRate0*theta/L) -
 	a*(1.0 - slipRate/slipRateLinear);
     } // else
-
     friction = -mu_f * normalTraction + properties[p_cohesion];
+    
+  } else {
+    friction = properties[p_cohesion];
+  } // if/else
 
-#if 0
-    std::cout << "slip: " << slip
-	      << ", slipRate: " << slipRate
-	      << ", stateVar: " << theta
-	      << ", mu_f: " << mu_f
-	      << std::endl;
-#endif
-
-#endif
-  } // if
-
   PetscLogFlops(12);
 
   return friction;
 } // _calcFriction
 
-// ----------------------------------------------------------------------
-// Compute change in friction for a change in slip (Jacobian).
-PylithScalar
-pylith::friction::RateStateAgeing::_calcFrictionSlope(const PylithScalar slip,
-						      const PylithScalar slipRate,
-						      const PylithScalar normalTraction,
-						      const PylithScalar* properties,
-						      const int numProperties,
-						      const PylithScalar* stateVars,
-						      const int numStateVars)
-{ // _calcFrictionSlope
-  assert(properties);
-  assert(_RateStateAgeing::numProperties == numProperties);
-  assert(numStateVars);
-  assert(_RateStateAgeing::numStateVars == numStateVars);
 
-  PylithScalar slope = 0.0;
-  if (normalTraction <= 0.0) {
-    // if fault is in compression
-
-    const PylithScalar a = properties[p_a];
-    const PylithScalar slipRate0 = properties[p_slipRate0];
-
-    const PylithScalar slipRateLinear = _minSlipRate;
-    const PylithScalar slipRateEff = std::max(slipRate, slipRateLinear);
-
-    //slope = -slipRateEff / (slipRate0 * normalTraction * a);
-    slope = -normalTraction * a; // log space
-  } // if
-
-  PetscLogFlops(5);
-
-  return slope;
-} // _calcFrictionSlope
-
-
 // ----------------------------------------------------------------------
 // Update state variables (for next time step).
 void

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.hh	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.hh	2012-02-24 01:10:59 UTC (rev 19668)
@@ -143,26 +143,6 @@
 			     const PylithScalar* stateVars,
 			     const int numStateVars);
 
-  /** Compute change in friction for a change in slip (Jacobian).
-   *
-   * @param slip Current slip at location.
-   * @param slipRate Current slip rate at location.
-   * @param normalTraction Normal traction at location.
-   * @param properties Properties at location.
-   * @param numProperties Number of properties.
-   * @param stateVars State variables at location.
-   * @param numStateVars Number of state variables.
-   *
-   * @returns Change in friction for a chance in slip (dT/dD).
-   */
-  PylithScalar _calcFrictionSlope(const PylithScalar slip,
-			    const PylithScalar slipRate,
-			    const PylithScalar normalTraction,
-			    const PylithScalar* properties,
-			    const int numProperties,
-			    const PylithScalar* stateVars,
-			    const int numStateVars);
-  
   /** Update state variables (for next time step).
    *
    * @param t Time in simulation.

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "pylith/materials/Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -305,43 +305,6 @@
 } // _calcFriction
 
 // ----------------------------------------------------------------------
-// Compute change in friction for a change in slip (Jacobian).
-PylithScalar
-pylith::friction::SlipWeakening::_calcFrictionSlope(const PylithScalar slip,
-						    const PylithScalar slipRate,
-						    const PylithScalar normalTraction,
-						    const PylithScalar* properties,
-						    const int numProperties,
-						    const PylithScalar* stateVars,
-						    const int numStateVars)
-{ // _calcFrictionSlope
-  assert(properties);
-  assert(_SlipWeakening::numProperties == numProperties);
-  assert(stateVars);
-  assert(_SlipWeakening::numStateVars == numStateVars);
-
-  PylithScalar slope = 0.0;
-  if (normalTraction <= 0.0) {
-    // if fault is in compression
-    const PylithScalar slipPrev = stateVars[s_slipPrev];
-    const PylithScalar slipCum = stateVars[s_slipCum] + fabs(slip - slipPrev);
-
-    if (slipCum < properties[p_d0]) {
-      // if/else linear slip-weakening form of mu_f 
-      slope = -normalTraction * (properties[p_coefS] - properties[p_coefD]) 
-	/ properties[p_d0];
-      } else {
-      slope = 0.0;
-      } // if/else
-  } // if
-
-  PetscLogFlops(7);
-
-  return slope;
-} // _calcFrictionSlope
-
-
-// ----------------------------------------------------------------------
 // Update state variables (for next time step).
 void
 pylith::friction::SlipWeakening::_updateStateVars(const PylithScalar t,

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh	2012-02-24 01:10:59 UTC (rev 19668)
@@ -131,26 +131,6 @@
 			     const PylithScalar* stateVars,
 			     const int numStateVars);
   
-  /** Compute change in friction for a change in slip (Jacobian).
-   *
-   * @param slip Current slip at location.
-   * @param slipRate Current slip rate at location.
-   * @param normalTraction Normal traction at location.
-   * @param properties Properties at location.
-   * @param numProperties Number of properties.
-   * @param stateVars State variables at location.
-   * @param numStateVars Number of state variables.
-   *
-   * @returns Change in friction for a chance in slip (dT/dD).
-   */
-  PylithScalar _calcFrictionSlope(const PylithScalar slip,
-				  const PylithScalar slipRate,
-				  const PylithScalar normalTraction,
-				  const PylithScalar* properties,
-				  const int numProperties,
-				  const PylithScalar* stateVars,
-				  const int numStateVars);
-  
   /** Update state variables (for next time step).
    *
    * @param t Time in simulation.

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "pylith/materials/Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -312,19 +312,6 @@
 } // _calcFriction
 
 // ----------------------------------------------------------------------
-// Compute change in friction for a change in slip (Jacobian).
-PylithScalar
-pylith::friction::SlipWeakeningTime::_calcFrictionSlope(const PylithScalar slip,
-							const PylithScalar slipRate,
-							const PylithScalar normalTraction,
-							const PylithScalar* properties,
-							const int numProperties,
-							const PylithScalar* stateVars,
-							const int numStateVars)
-{ // _calcFrictionSlope
-} // _calcFrictionSlope
-
-// ----------------------------------------------------------------------
 // Update state variables (for next time step).
 void
 pylith::friction::SlipWeakeningTime::_updateStateVars(const PylithScalar t,

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.hh	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.hh	2012-02-24 01:10:59 UTC (rev 19668)
@@ -125,26 +125,6 @@
 			     const PylithScalar* stateVars,
 			     const int numStateVars);
 
-  /** Compute change in friction for a change in slip (Jacobian).
-   *
-   * @param slip Current slip at location.
-   * @param slipRate Current slip rate at location.
-   * @param normalTraction Normal traction at location.
-   * @param properties Properties at location.
-   * @param numProperties Number of properties.
-   * @param stateVars State variables at location.
-   * @param numStateVars Number of state variables.
-   *
-   * @returns Change in friction for a chance in slip (dT/dD).
-   */
-  PylithScalar _calcFrictionSlope(const PylithScalar slip,
-				  const PylithScalar slipRate,
-				  const PylithScalar normalTraction,
-				  const PylithScalar* properties,
-				  const int numProperties,
-				  const PylithScalar* stateVars,
-				  const int numStateVars);
-  
   /** Update state variables (for next time step).
    *
    * @param t Time in simulation.

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "pylith/materials/Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -159,7 +159,7 @@
 
   const PylithScalar friction = (normalTraction <= 0.0) ?
     -properties[p_coef] * normalTraction + properties[p_cohesion]: 
-    0.0;
+    properties[p_cohesion];
 
   PetscLogFlops(2);
 
@@ -167,21 +167,4 @@
 } // _calcFriction
 
 
-// ----------------------------------------------------------------------
-// Compute change in friction for a change in slip (Jacobian).
-PylithScalar
-pylith::friction::StaticFriction::_calcFrictionSlope(const PylithScalar slip,
-						     const PylithScalar slipRate,
-						     const PylithScalar normalTraction,
-						     const PylithScalar* properties,
-						     const int numProperties,
-						     const PylithScalar* stateVars,
-						     const int numStateVars)
-{ // _calcFrictionSlope
-  const PylithScalar slope = 0;
-
-  return slope;
-} // _calcFrictionSlope
-
-
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.hh	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.hh	2012-02-24 01:10:59 UTC (rev 19668)
@@ -96,26 +96,6 @@
 			     const PylithScalar* stateVars,
 			     const int numStateVars);
 
-  /** Compute change in friction for a change in slip (Jacobian).
-   *
-   * @param slip Current slip at location.
-   * @param slipRate Current slip rate at location.
-   * @param normalTraction Normal traction at location.
-   * @param properties Properties at location.
-   * @param numProperties Number of properties.
-   * @param stateVars State variables at location.
-   * @param numStateVars Number of state variables.
-   *
-   * @returns Change in friction for a chance in slip (dT/dD).
-   */
-  PylithScalar _calcFrictionSlope(const PylithScalar slip,
-			    const PylithScalar slipRate,
-			    const PylithScalar normalTraction,
-			    const PylithScalar* properties,
-			    const int numProperties,
-			    const PylithScalar* stateVars,
-			    const int numStateVars);
-  
   // PRIVATE MEMBERS ////////////////////////////////////////////////////
 private :
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "pylith/materials/Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -274,7 +274,9 @@
 	mu_f = properties[p_coefD];
       } // if/else
     friction = - mu_f * normalTraction + properties[p_cohesion];
-  } // if
+  } else {
+    friction = properties[p_cohesion];
+  } // if/else
 
   PetscLogFlops(6);
 
@@ -310,21 +312,4 @@
 } // _updateStateVars
 
 
-// ----------------------------------------------------------------------
-// Compute change in friction for a change in slip (Jacobian).
-PylithScalar
-pylith::friction::TimeWeakening::_calcFrictionSlope(const PylithScalar slip,
-						    const PylithScalar slipRate,
-						    const PylithScalar normalTraction,
-						    const PylithScalar* properties,
-						    const int numProperties,
-						    const PylithScalar* stateVars,
-						    const int numStateVars)
-{ // _calcFrictionSlope
-  const PylithScalar slope = 0;
-
-  return slope;
-} // _calcFrictionSlope
-
-
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.hh	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.hh	2012-02-24 01:10:59 UTC (rev 19668)
@@ -125,26 +125,6 @@
 			     const PylithScalar* stateVars,
 			     const int numStateVars);
   
-  /** Compute change in friction for a change in slip (Jacobian).
-   *
-   * @param slip Current slip at location.
-   * @param slipRate Current slip rate at location.
-   * @param normalTraction Normal traction at location.
-   * @param properties Properties at location.
-   * @param numProperties Number of properties.
-   * @param stateVars State variables at location.
-   * @param numStateVars Number of state variables.
-   *
-   * @returns Change in friction for a chance in slip (dT/dD).
-   */
-  PylithScalar _calcFrictionSlope(const PylithScalar slip,
-			    const PylithScalar slipRate,
-			    const PylithScalar normalTraction,
-			    const PylithScalar* properties,
-			    const int numProperties,
-			    const PylithScalar* stateVars,
-			    const int numStateVars);
-  
   /** Update state variables (for next time step).
    *
    * @param slip Current slip at location.

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -415,7 +415,7 @@
 
   // It's unclear what to do for an elasto-plastic material, which has no
   // inherent time scale. For now, just set dtStable to a large value.
-  const PylithScalar dtStable = pylith::PYLITH_MAXDOUBLE;
+  const PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
 
   return dtStable;
 } // _stableTimeStepImplicit

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -435,7 +435,7 @@
   assert(_numVarsQuadPt == numStateVars);
   // It's unclear what to do for an elasto-plastic material, which has no
   // inherent time scale. For now, just set dtStable to a large value.
-  const PylithScalar dtStable = pylith::PYLITH_MAXDOUBLE;
+  const PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
 
   return dtStable;
 } // _stableTimeStepImplicit

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -327,7 +327,7 @@
 PylithScalar
 pylith::materials::ElasticIsotropic3D::stableTimeStepImplicit(
 					const topology::Mesh& mesh) {
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 }
 
 // ----------------------------------------------------------------------
@@ -339,7 +339,7 @@
 				     const PylithScalar* stateVars,
 				     const int numStateVars) const
 { // _stableTimeStepImplicit
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 } // _stableTimeStepImplicit
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "pylith/topology/Mesh.hh" // USES Mesh
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/utils/array.hh" // USES scalar_array, std::vector
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -243,7 +243,7 @@
   assert(_initialStressCell.size() == numQuadPts*_tensorSize);
   assert(_initialStrainCell.size() == numQuadPts*_tensorSize);
 
-  PylithScalar dtStable = pylith::PYLITH_MAXDOUBLE;
+  PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
 
   // Get cells associated with material
   const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticPlaneStrain.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticPlaneStrain.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -313,7 +313,7 @@
 PylithScalar
 pylith::materials::ElasticPlaneStrain::stableTimeStepImplicit(
 					const topology::Mesh& mesh) {
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 }
 
 // ----------------------------------------------------------------------
@@ -325,7 +325,7 @@
 				     const PylithScalar* stateVars,
 				     const int numStateVars) const
 { // _stableTimeStepImplicit
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 } // _stableTimeStepImplicit
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticPlaneStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticPlaneStress.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticPlaneStress.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -316,7 +316,7 @@
 PylithScalar
 pylith::materials::ElasticPlaneStress::stableTimeStepImplicit(
 					const topology::Mesh& mesh) {
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 }
 
 // ----------------------------------------------------------------------
@@ -328,7 +328,7 @@
 				     const PylithScalar* stateVars,
 				     const int numStateVars) const
 { // _stableTimeStepImplicit
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 } // _stableTimeStepImplicit
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticStrain1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticStrain1D.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticStrain1D.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -293,7 +293,7 @@
 PylithScalar
 pylith::materials::ElasticStrain1D::stableTimeStepImplicit(
 					const topology::Mesh& mesh) {
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 }
 
 // ----------------------------------------------------------------------
@@ -305,7 +305,7 @@
 					      const PylithScalar* stateVars,
 					      const int numStateVars) const
 { // _stableTimeStepImplicit
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 } // _stableTimeStepImplicit
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticStress1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticStress1D.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticStress1D.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,7 +23,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -293,7 +293,7 @@
 PylithScalar
 pylith::materials::ElasticStress1D::stableTimeStepImplicit(
 					const topology::Mesh& mesh) {
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 }
 
 // ----------------------------------------------------------------------
@@ -305,7 +305,7 @@
 					      const PylithScalar* stateVars,
 					      const int numStateVars) const
 { // _stableTimeStepImplicit
-  return pylith::PYLITH_MAXDOUBLE;
+  return pylith::PYLITH_MAXSCALAR;
 } // _stableTimeStepImplicit
 
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -24,7 +24,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -294,7 +294,7 @@
     PylithScalar muRatio = dbValues[db_shearRatio + imodel];
     PylithScalar viscosity = dbValues[db_viscosity + imodel];
     PylithScalar muFac = muRatio*mu;
-    PylithScalar maxwellTime = pylith::PYLITH_MAXDOUBLE;
+    PylithScalar maxwellTime = pylith::PYLITH_MAXSCALAR;
     if (muFac > 0.0)
       maxwellTime = viscosity / muFac;
     if (muRatio < 0.0 || viscosity < 0.0 || muFac < 0.0 || maxwellTime < 0.0) {
@@ -697,7 +697,7 @@
   PylithScalar shearRatio = 0.0;
   for (int imodel = 0; imodel < numMaxwellModels; ++imodel) {
     shearRatio = properties[p_shearRatio + imodel];
-    PylithScalar maxwellTime = pylith::PYLITH_MAXDOUBLE;
+    PylithScalar maxwellTime = pylith::PYLITH_MAXSCALAR;
     visFrac += shearRatio;
     if (shearRatio != 0.0) {
       maxwellTime = properties[p_maxwellTime + imodel];
@@ -882,7 +882,7 @@
   assert(0 != stateVars);
   assert(_numVarsQuadPt == numStateVars);
 
-  PylithScalar dtStable = pylith::PYLITH_MAXDOUBLE;
+  PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
 
   const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
   for (int i=0; i < numMaxwellModels; ++i) {

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -24,7 +24,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -289,7 +289,7 @@
     PylithScalar muRatio = dbValues[db_shearRatio + imodel];
     PylithScalar viscosity = dbValues[db_viscosity + imodel];
     PylithScalar muFac = muRatio * mu;
-    PylithScalar maxwellTime = pylith::PYLITH_MAXDOUBLE;
+    PylithScalar maxwellTime = pylith::PYLITH_MAXSCALAR;
     if (muFac > 0.0)
       maxwellTime = viscosity / muFac;
     if (muRatio < 0.0 || viscosity < 0.0 || muFac < 0.0 || maxwellTime < 0.0) {
@@ -685,7 +685,7 @@
   PylithScalar shearRatio = 0.0;
   for (int imodel = 0; imodel < numMaxwellModels; ++imodel) {
     shearRatio = properties[p_shearRatio + imodel];
-    PylithScalar maxwellTime = pylith::PYLITH_MAXDOUBLE;
+    PylithScalar maxwellTime = pylith::PYLITH_MAXSCALAR;
     visFrac += shearRatio;
     if (shearRatio != 0.0) {
       maxwellTime = properties[p_maxwellTime + imodel];
@@ -830,7 +830,7 @@
   assert(0 != stateVars);
   assert(_numVarsQuadPt == numStateVars);
 
-  PylithScalar dtStable = pylith::PYLITH_MAXDOUBLE;
+  PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
 
   const int numMaxwellModels = _GenMaxwellPlaneStrain::numMaxwellModels;
   for (int i=0; i < numMaxwellModels; ++i) {

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellQpQsIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellQpQsIsotropic3D.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellQpQsIsotropic3D.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -24,7 +24,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -350,8 +350,8 @@
     PylithScalar bulkRatio = dbValues[db_bulkRatio + imodel];
     PylithScalar shearViscosity = dbValues[db_shearViscosity + imodel];
     PylithScalar bulkViscosity = dbValues[db_bulkViscosity + imodel];
-    PylithScalar maxwellTimeShear = pylith::PYLITH_MAXDOUBLE;
-    PylithScalar maxwellTimeBulk = pylith::PYLITH_MAXDOUBLE;
+    PylithScalar maxwellTimeShear = pylith::PYLITH_MAXSCALAR;
+    PylithScalar maxwellTimeBulk = pylith::PYLITH_MAXSCALAR;
     maxwellTimeShear = shearViscosity / mu;
     maxwellTimeBulk = bulkViscosity / k;
     if (shearRatio < 0.0 || shearViscosity < 0.0 || maxwellTimeShear < 0.0 || 
@@ -954,7 +954,7 @@
 
   const int numMaxwellModels = _GenMaxwellQpQsIsotropic3D::numMaxwellModels;
 
-  PylithScalar dtStable = pylith::PYLITH_MAXDOUBLE;
+  PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
 
   for (int i=0; i < numMaxwellModels; ++i) {
     const PylithScalar maxwellTime = properties[p_maxwellTimeShear+i];

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/MaxwellPlaneStrain.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -24,7 +24,7 @@
 #include "Metadata.hh" // USES Metadata
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/PowerLaw3D.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/PowerLaw3D.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -24,7 +24,7 @@
 #include "EffectiveStress.hh" // USES EffectiveStress
 
 #include "pylith/utils/array.hh" // USES scalar_array
-#include "pylith/utils/constdefs.h" // USES PYLITH_MAXDOUBLE
+#include "pylith/utils/constdefs.h" // USES PYLITH_MAXSCALAR
 
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
@@ -420,7 +420,7 @@
   const PylithScalar effStress = sqrt(0.5 * devStressProd);
   PylithScalar dtTest = 0.0;
   if (effStress <= 0.0) {
-    dtTest = pylith::PYLITH_MAXDOUBLE;
+    dtTest = pylith::PYLITH_MAXSCALAR;
   } else {
     dtTest = 0.05 *
     pow((referenceStress/effStress), (powerLawExp - 1.0)) *

Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,6 +23,7 @@
 #include "Formulation.hh" // USES Formulation
 #include "pylith/topology/SolutionFields.hh" // USES SolutionFields
 #include "pylith/topology/Jacobian.hh" // USES Jacobian
+#include "pylith/utils/constdefs.h" // USES PYLITH_MAXSCALAR
 #include "pylith/utils/EventLogger.hh" // USES EventLogger
 
 #include <petscsnes.h> // USES PetscSNES
@@ -243,7 +244,10 @@
     *ynorm = snes->maxstep;
   }
   ierr      = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr);
-  minlambda = snes->steptol/rellength;
+
+  // Place reasonable absolute limit on minimum lambda
+  minlambda = std::max(snes->steptol/rellength, 1.0/PYLITH_MAXSCALAR);
+
   ierr      = MatMult(snes->jacobian,y,w);CHKERRQ(ierr);
 #if defined(PETSC_USE_COMPLEX)
   ierr      = VecDot(f,w,&cinitslope);CHKERRQ(ierr);
@@ -325,10 +329,12 @@
 	ierr = PetscViewerASCIIPrintf(snes->ls_monitor,"    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",(double)fnorm,(double)*gnorm,(double)*ynorm,(double)minlambda,(double)lambda,(double)initslope);CHKERRQ(ierr);
         ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
       }
+#if 0 // DEBUGGING
       assert(lsctx);
       Formulation* formulation = (Formulation*) lsctx;
       assert(formulation);
       formulation->printState(&w, &g, &x, &y);
+#endif
 #if 0 // Original PETSc code
       *flag = PETSC_FALSE; // DIVERGED_LINE_SEARCH
 #else

Modified: short/3D/PyLith/trunk/libsrc/pylith/utils/constdefs.h
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/utils/constdefs.h	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/libsrc/pylith/utils/constdefs.h	2012-02-24 01:10:59 UTC (rev 19668)
@@ -28,9 +28,9 @@
 #include "types.hh" // HASA PylithScalar
 
 namespace pylith {
-  static const double PYLITH_MAXDOUBLE = 1.0e+30;
+  static const double PYLITH_MAXDOUBLE = 1.0e+99;
   static const float PYLITH_MAXFLOAT = 1.0e+30;
-  static const PylithScalar PYLITH_MAXSCALAR = 1.0e+30;
+  static const PylithScalar PYLITH_MAXSCALAR = (sizeof(PylithScalar) == sizeof(double)) ? PYLITH_MAXDOUBLE : PYLITH_MAXSCALAR;
 }
     
 

Modified: short/3D/PyLith/trunk/modulesrc/friction/FrictionModel.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/friction/FrictionModel.i	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/modulesrc/friction/FrictionModel.i	2012-02-24 01:10:59 UTC (rev 19668)
@@ -149,21 +149,6 @@
 				const PylithScalar normalTraction);
   
 
-      /** Compute change in friction for a change in slip (Jacobian).
-       *
-       * @pre Must call retrievePropsAndVars for cell before calling
-       * calcFriction().
-       *
-       * @param slip Current slip at location.
-       * @param slipRate Current slip rate at location.
-       * @param normalTraction Normal traction at location.
-       *
-       * @returns Change in friction for a chance in slip (dT/dD).
-       */
-      PylithScalar calcFrictionSlope(const PylithScalar slip,
-			       const PylithScalar slipRate,
-			       const PylithScalar normalTraction);
-  
       /** Compute friction at vertex.
        *
        * @pre Must call retrievePropsAndVars for cell before calling
@@ -261,27 +246,6 @@
 				 const PylithScalar* stateVars,
 				 const int numStateVars) = 0;
 
-      /** Compute change in friction for a change in slip (Jacobian).
-       *
-       * @param slip Current slip at location.
-       * @param slipRate Current slip rate at location.
-       * @param normalTraction Normal traction at location.
-       * @param properties Properties at location.
-       * @param numProperties Number of properties.
-       * @param stateVars State variables at location.
-       * @param numStateVars Number of state variables.
-       *
-       * @returns Change in friction for a chance in slip (dT/dD).
-       */
-      virtual
-      PylithScalar _calcFrictionSlope(const PylithScalar slip,
-				const PylithScalar slipRate,
-				const PylithScalar normalTraction,
-				const PylithScalar* properties,
-				const int numProperties,
-				const PylithScalar* stateVars,
-				const int numStateVars) = 0;
-  
       /** Update state variables (for next time step).
        *
        * @param t Current time.

Modified: short/3D/PyLith/trunk/modulesrc/friction/RateStateAgeing.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/friction/RateStateAgeing.i	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/modulesrc/friction/RateStateAgeing.i	2012-02-24 01:10:59 UTC (rev 19668)
@@ -90,26 +90,6 @@
 				 const PylithScalar* stateVars,
 				 const int numStateVars);
 
-      /** Compute change in friction for a change in slip (Jacobian).
-       *
-       * @param slip Current slip at location.
-       * @param slipRate Current slip rate at location.
-       * @param normalTraction Normal traction at location.
-       * @param properties Properties at location.
-       * @param numProperties Number of properties.
-       * @param stateVars State variables at location.
-       * @param numStateVars Number of state variables.
-       *
-       * @returns Change in friction for a chance in slip (dT/dD).
-       */
-      PylithScalar _calcFrictionSlope(const PylithScalar slip,
-				const PylithScalar slipRate,
-				const PylithScalar normalTraction,
-				const PylithScalar* properties,
-				const int numProperties,
-				const PylithScalar* stateVars,
-				const int numStateVars);
-
       /** Update state variables (for next time step).
        *
        * @param slip Current slip at location.

Modified: short/3D/PyLith/trunk/modulesrc/friction/SlipWeakening.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/friction/SlipWeakening.i	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/modulesrc/friction/SlipWeakening.i	2012-02-24 01:10:59 UTC (rev 19668)
@@ -89,26 +89,6 @@
 				 const PylithScalar* stateVars,
 				 const int numStateVars);
 
-      /** Compute change in friction for a change in slip (Jacobian).
-       *
-       * @param slip Current slip at location.
-       * @param slipRate Current slip rate at location.
-       * @param normalTraction Normal traction at location.
-       * @param properties Properties at location.
-       * @param numProperties Number of properties.
-       * @param stateVars State variables at location.
-       * @param numStateVars Number of state variables.
-       *
-       * @returns Change in friction for a chance in slip (dT/dD).
-       */
-      PylithScalar _calcFrictionSlope(const PylithScalar slip,
-				const PylithScalar slipRate,
-				const PylithScalar normalTraction,
-				const PylithScalar* properties,
-				const int numProperties,
-				const PylithScalar* stateVars,
-				const int numStateVars);
-  
       /** Update state variables (for next time step).
        *
        * @param slip Current slip at location.

Modified: short/3D/PyLith/trunk/modulesrc/friction/SlipWeakeningTime.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/friction/SlipWeakeningTime.i	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/modulesrc/friction/SlipWeakeningTime.i	2012-02-24 01:10:59 UTC (rev 19668)
@@ -83,26 +83,6 @@
 				 const PylithScalar* stateVars,
 				 const int numStateVars);
 
-      /** Compute change in friction for a change in slip (Jacobian).
-       *
-       * @param slip Current slip at location.
-       * @param slipRate Current slip rate at location.
-       * @param normalTraction Normal traction at location.
-       * @param properties Properties at location.
-       * @param numProperties Number of properties.
-       * @param stateVars State variables at location.
-       * @param numStateVars Number of state variables.
-       *
-       * @returns Change in friction for a chance in slip (dT/dD).
-       */
-      PylithScalar _calcFrictionSlope(const PylithScalar slip,
-				const PylithScalar slipRate,
-				const PylithScalar normalTraction,
-				const PylithScalar* properties,
-				const int numProperties,
-				const PylithScalar* stateVars,
-				const int numStateVars);
-  
     }; // class SlipWeakeningTime
 
   } // friction

Modified: short/3D/PyLith/trunk/modulesrc/friction/StaticFriction.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/friction/StaticFriction.i	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/modulesrc/friction/StaticFriction.i	2012-02-24 01:10:59 UTC (rev 19668)
@@ -83,26 +83,6 @@
 				 const PylithScalar* stateVars,
 				 const int numStateVars);
 
-      /** Compute change in friction for a change in slip (Jacobian).
-       *
-       * @param slip Current slip at location.
-       * @param slipRate Current slip rate at location.
-       * @param normalTraction Normal traction at location.
-       * @param properties Properties at location.
-       * @param numProperties Number of properties.
-       * @param stateVars State variables at location.
-       * @param numStateVars Number of state variables.
-       *
-       * @returns Change in friction for a chance in slip (dT/dD).
-       */
-      PylithScalar _calcFrictionSlope(const PylithScalar slip,
-				const PylithScalar slipRate,
-				const PylithScalar normalTraction,
-				const PylithScalar* properties,
-				const int numProperties,
-				const PylithScalar* stateVars,
-				const int numStateVars);
-  
     }; // class StaticFriction
 
   } // friction

Modified: short/3D/PyLith/trunk/modulesrc/friction/TimeWeakening.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/friction/TimeWeakening.i	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/modulesrc/friction/TimeWeakening.i	2012-02-24 01:10:59 UTC (rev 19668)
@@ -83,26 +83,6 @@
 				 const PylithScalar* stateVars,
 				 const int numStateVars);
 
-      /** Compute change in friction for a change in slip (Jacobian).
-       *
-       * @param slip Current slip at location.
-       * @param slipRate Current slip rate at location.
-       * @param normalTraction Normal traction at location.
-       * @param properties Properties at location.
-       * @param numProperties Number of properties.
-       * @param stateVars State variables at location.
-       * @param numStateVars Number of state variables.
-       *
-       * @returns Change in friction for a chance in slip (dT/dD).
-       */
-      PylithScalar _calcFrictionSlope(const PylithScalar slip,
-				const PylithScalar slipRate,
-				const PylithScalar normalTraction,
-				const PylithScalar* properties,
-				const int numProperties,
-				const PylithScalar* stateVars,
-				const int numStateVars);
-  
       /** Update state variables (for next time step).
        *
        * @param slip Current slip at location.

Modified: short/3D/PyLith/trunk/modulesrc/utils/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/utils/Makefile.am	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/modulesrc/utils/Makefile.am	2012-02-24 01:10:59 UTC (rev 19668)
@@ -41,7 +41,8 @@
 	utils.i \
 	pylith_general.i \
 	EventLogger.i \
-	TestArray.i
+	TestArray.i \
+	constdefs.i
 
 utils_swig_generated = \
 	utils_wrap.cxx \

Copied: short/3D/PyLith/trunk/modulesrc/utils/constdefs.i (from rev 19667, short/3D/PyLith/branches/v1.6-stable/modulesrc/utils/constdefs.i)
===================================================================
--- short/3D/PyLith/trunk/modulesrc/utils/constdefs.i	                        (rev 0)
+++ short/3D/PyLith/trunk/modulesrc/utils/constdefs.i	2012-02-24 01:10:59 UTC (rev 19668)
@@ -0,0 +1,58 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2011 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file modulesrc/utils/constdefs.i
+ *
+ * @brief PyLith constants.
+ */
+
+// ----------------------------------------------------------------------
+// PYLITH_MAXDOUBLE
+%inline %{
+  double
+  maxdouble(void)
+  { // maxdouble
+    return pylith::PYLITH_MAXDOUBLE;
+  } // maxdouble
+%} // inline
+
+
+// ----------------------------------------------------------------------
+// PYLITH_MAXFLOAT
+%inline %{
+  float
+  maxfloat(void)
+  { // maxfloat
+    return pylith::PYLITH_MAXFLOAT;
+  } // maxfloat
+%} // inline
+
+
+// ----------------------------------------------------------------------
+// PYLITH_MAXSCALAR
+%inline %{
+  double
+  maxscalar(void)
+  { // maxscalar
+    return pylith::PYLITH_MAXSCALAR;
+  } // maxscalar
+%} // inline
+
+
+// End of file 

Modified: short/3D/PyLith/trunk/modulesrc/utils/utils.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/utils/utils.i	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/modulesrc/utils/utils.i	2012-02-24 01:10:59 UTC (rev 19668)
@@ -23,6 +23,7 @@
 %{
 #include "pylith/utils/EventLogger.hh"
 #include "pylith/utils/TestArray.hh"
+#include "pylith/utils/constdefs.h"
 
 #include <petsclog.h> // USES PetscLogEventBegin/End() in inline methods
 #include "pylith/utils/arrayfwd.hh" // USES scalar_array
@@ -56,6 +57,7 @@
 %include "pylith_general.i"
 %include "EventLogger.i"
 %include "TestArray.i"
+%include "constdefs.i"
 
 // End of file
 

Modified: short/3D/PyLith/trunk/tests/2d/frictionslide/plot_friction.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/plot_friction.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/plot_friction.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 
-sim = "ratestate_weak"
+sim = "ratestate_stable"
 
 # ======================================================================
 import tables

Modified: short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg	2012-02-24 01:10:59 UTC (rev 19668)
@@ -78,7 +78,7 @@
 
 [pylithapp.timedependent.interfaces]
 fault = pylith.faults.FaultCohesiveDyn
-fault.zero_tolerance = 1.0e-12
+fault.zero_tolerance = 1.0e-13
 
 [pylithapp.timedependent.interfaces.fault]
 id = 100

Modified: short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg	2012-02-24 01:10:59 UTC (rev 19668)
@@ -33,7 +33,7 @@
 # problem
 # ----------------------------------------------------------------------
 [pylithapp.timedependent.formulation.time_step]
-total_time = 60.0*hour ; total time of simulation
+total_time = 90.0*hour ; total time of simulation
 dt = 1.0*hour
 
 [pylithapp.timedependent.normalizer]
@@ -74,7 +74,7 @@
 db_change = spatialdata.spatialdb.UniformDB
 db_change.label = Amplitude of Dirichlet BC on +x
 db_change.values = [displacement-x,displacement-y,displacement-z,change-start-time]
-db_change.data = [-0.5*m,-0.2*m,0.0*m,0.0*hour]
+db_change.data = [-0.5*m,-0.3*m,0.0*m,0.0*hour]
 
 th_change = spatialdata.spatialdb.TimeHistory
 th_change.label = Time history of tidal load
@@ -89,7 +89,7 @@
 db_change = spatialdata.spatialdb.UniformDB
 db_change.label = Amplitude of Dirichlet BC on -x
 db_change.values = [displacement-x,displacement-y,displacement-z,change-start-time]
-db_change.data = [+0.5*m,+0.2*m,0.0*m,0.0*hour]
+db_change.data = [+0.5*m,+0.3*m,0.0*m,0.0*hour]
 
 th_change = spatialdata.spatialdb.TimeHistory
 th_change.label = Time history of tidal load
@@ -114,7 +114,7 @@
 
 [pylithapp.timedependent.interfaces.fault]
 label = fault
-zero_tolerance = 1.0e-10
+zero_tolerance = 1.0e-12
 
 friction = pylith.friction.SlipWeakening
 friction.label = Slip weakening
@@ -171,7 +171,7 @@
 sub_pc_factor_shift_type = nonzero
 
 # Convergence parameters.
-ksp_rtol = 1.0e-10
+ksp_rtol = 1.0e-12
 ksp_atol = 1.0e-12
 ksp_max_it = 100
 ksp_gmres_restart = 50
@@ -182,11 +182,12 @@
 ksp_converged_reason = true
 
 # Nonlinear solver monitoring options.
-snes_rtol = 1.0e-10
-snes_atol = 1.0e-9
+snes_rtol = 1.0e-12
+snes_atol = 1.0e-10
 snes_max_it = 100
 snes_monitor = true
 snes_ls_monitor = true
+#snes_steptol = 1.0e-20
 #snes_view = true
 snes_converged_reason = true
 snes_monitor_solution_update = true

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -1086,6 +1086,246 @@
   4.1,  5.2,  5.3,
   5.4,  5.5,  5.6,
   7.7,  5.8,  5.9,
+  0.0,  0.0,  0.0, // 18x
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 18y
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 18z
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 19x
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 19y
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 19z
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 20x
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 20y
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 20z
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 21x
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 21y
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0, // 21z
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
+  0.0,  0.0,  0.0,
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataHex8::_orientation[] = {

Modified: short/3D/PyLith/trunk/unittests/libtests/friction/data/StaticFrictionData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/friction/data/StaticFrictionData.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/friction/data/StaticFrictionData.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -85,7 +85,7 @@
 
 const PylithScalar pylith::friction::StaticFrictionData::_friction[] = {
   1000001.32,
-  0.0,
+  1.0e+6,
 };
 
 const PylithScalar pylith::friction::StaticFrictionData::_slip[] = {

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDep.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -189,7 +189,7 @@
                                                plasStrainB,
                                                initialStressB, initialStrainB)
 
-    self.dtStableImplicit = 1.0e30
+    self.dtStableImplicit = 1.0e+99
 
     return
 

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPrager3DTimeDepData.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -45,7 +45,7 @@
 
 const PylithScalar pylith::materials::DruckerPrager3DTimeDepData::_densityScale =   1.00000000e+03;
 
-const PylithScalar pylith::materials::DruckerPrager3DTimeDepData::_dtStableImplicit =   1.00000000e+30;
+const PylithScalar pylith::materials::DruckerPrager3DTimeDepData::_dtStableImplicit =   1.00000000e+99;
 
 const int pylith::materials::DruckerPrager3DTimeDepData::_numPropertyValues[] = {
 1,

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/DruckerPragerPlaneStrainTimeDepData.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -45,7 +45,7 @@
 
 const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_densityScale =   1.00000000e+03;
 
-const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_dtStableImplicit =   1.00000000e+30;
+const PylithScalar pylith::materials::DruckerPragerPlaneStrainTimeDepData::_dtStableImplicit =   1.00000000e+99;
 
 const int pylith::materials::DruckerPragerPlaneStrainTimeDepData::_numPropertyValues[] = {
 1,

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -45,7 +45,7 @@
 
 const PylithScalar pylith::materials::ElasticIsotropic3DData::_densityScale =   1.00000000e+03;
 
-const PylithScalar pylith::materials::ElasticIsotropic3DData::_dtStableImplicit =   1.00000000e+30;
+const PylithScalar pylith::materials::ElasticIsotropic3DData::_dtStableImplicit =   1.00000000e+99;
 
 const int pylith::materials::ElasticIsotropic3DData::_numPropertyValues[] = {
 1,

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -88,7 +88,7 @@
     self.densityScale = 0
 
     # Elastic material information
-    self.dtStableImplicit = 1.0e+30
+    self.dtStableImplicit = 1.0e+99
     self.density = None
     self.strain = None
     self.stress = None

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -45,7 +45,7 @@
 
 const PylithScalar pylith::materials::ElasticPlaneStrainData::_densityScale =   1.00000000e+03;
 
-const PylithScalar pylith::materials::ElasticPlaneStrainData::_dtStableImplicit =   1.00000000e+30;
+const PylithScalar pylith::materials::ElasticPlaneStrainData::_dtStableImplicit =   1.00000000e+99;
 
 const int pylith::materials::ElasticPlaneStrainData::_numPropertyValues[] = {
 1,

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -45,7 +45,7 @@
 
 const PylithScalar pylith::materials::ElasticPlaneStressData::_densityScale =   1.00000000e+03;
 
-const PylithScalar pylith::materials::ElasticPlaneStressData::_dtStableImplicit =   1.00000000e+30;
+const PylithScalar pylith::materials::ElasticPlaneStressData::_dtStableImplicit =   1.00000000e+99;
 
 const int pylith::materials::ElasticPlaneStressData::_numPropertyValues[] = {
 1,

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -45,7 +45,7 @@
 
 const PylithScalar pylith::materials::ElasticStrain1DData::_densityScale =   1.00000000e+03;
 
-const PylithScalar pylith::materials::ElasticStrain1DData::_dtStableImplicit =   1.00000000e+30;
+const PylithScalar pylith::materials::ElasticStrain1DData::_dtStableImplicit =   1.00000000e+99;
 
 const int pylith::materials::ElasticStrain1DData::_numPropertyValues[] = {
 1,

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc	2012-02-24 01:10:59 UTC (rev 19668)
@@ -45,7 +45,7 @@
 
 const PylithScalar pylith::materials::ElasticStress1DData::_densityScale =   1.00000000e+03;
 
-const PylithScalar pylith::materials::ElasticStress1DData::_dtStableImplicit =   1.00000000e+30;
+const PylithScalar pylith::materials::ElasticStress1DData::_dtStableImplicit =   1.00000000e+99;
 
 const int pylith::materials::ElasticStress1DData::_numPropertyValues[] = {
 1,

Modified: short/3D/PyLith/trunk/unittests/pytests/bc/TestAbsorbingDampers.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/TestAbsorbingDampers.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/TestAbsorbingDampers.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -80,7 +80,8 @@
     """
     (mesh, bc, fields) = self._initialize()
 
-    self.assertAlmostEqual(1.0, bc.stableTimeStep(mesh)/1.0e+30, 5)
+    from pylith.utils.utils import maxscalar
+    self.assertAlmostEqual(1.0, bc.stableTimeStep(mesh)/maxscalar(), 7)
     return
 
   

Modified: short/3D/PyLith/trunk/unittests/pytests/bc/TestNeumann.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/TestNeumann.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/TestNeumann.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -80,7 +80,8 @@
     """
     (mesh, bc, fields) = self._initialize()
 
-    self.assertAlmostEqual(1.0, bc.stableTimeStep(mesh)/1.0e+30, 5)
+    from pylith.utils.utils import maxscalar
+    self.assertAlmostEqual(1.0, bc.stableTimeStep(mesh)/maxscalar(), 7)
     return
 
   

Modified: short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveDyn.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveDyn.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -88,7 +88,7 @@
     fault.inventory.faultLabel = "fault"
     fault._configure()
 
-    nvertices = fault.numVertices(mesh)
+    nvertices = fault.numVerticesNoMesh(mesh)
     firstFaultVertex = 0
     firstLagrangeVertex = nvertices
     firstFaultCell      = 2*nvertices
@@ -132,7 +132,8 @@
     """
     (mesh, fault, fields) = self._initialize()
 
-    self.assertAlmostEqual(1.0, fault.stableTimeStep(mesh)/1.0e+30, 5)
+    from pylith.utils.utils import maxscalar
+    self.assertAlmostEqual(1.0, fault.stableTimeStep(mesh)/maxscalar(), 7)
     return
 
   
@@ -315,7 +316,7 @@
     fault.inventory.friction = friction
     fault._configure()
 
-    nvertices = fault.numVertices(mesh)
+    nvertices = fault.numVerticesNoMesh(mesh)
     firstFaultVertex = 0
     firstLagrangeVertex = nvertices
     firstFaultCell      = 2*nvertices

Modified: short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -117,7 +117,7 @@
     fault.inventory.faultLabel = "fault"
     fault._configure()
 
-    nvertices = fault.numVertices(mesh)
+    nvertices = fault.numVerticesNoMesh(mesh)
     firstFaultVertex = 0
     firstLagrangeVertex = nvertices
     firstFaultCell      = 2*nvertices
@@ -161,7 +161,8 @@
     """
     (mesh, fault, fields) = self._initialize()
 
-    self.assertAlmostEqual(1.0, fault.stableTimeStep(mesh)/1.0e+30, 5)
+    from pylith.utils.utils import maxscalar
+    self.assertAlmostEqual(1.0, fault.stableTimeStep(mesh)/maxscalar(), 7)
     return
 
   
@@ -347,7 +348,7 @@
     eqsrc.inventory.slipfn = slipfn
     eqsrc._configure()
 
-    nvertices = fault.numVertices(mesh)
+    nvertices = fault.numVerticesNoMesh(mesh)
     firstFaultVertex = 0
     firstLagrangeVertex = nvertices
     firstFaultCell      = 2*nvertices

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -121,7 +121,8 @@
     (mesh, integrator) = self._preinitialize()
     fields = self._initialize(mesh, integrator)
 
-    self.assertAlmostEqual(1.0, integrator.stableTimeStep(mesh)/1.0e+30, places=5)
+    from pylith.utils.utils import maxscalar
+    self.assertAlmostEqual(1.0, integrator.stableTimeStep(mesh)/maxscalar(), 7)
     return
 
 

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicitLgDeform.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicitLgDeform.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicitLgDeform.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -121,7 +121,8 @@
     (mesh, integrator) = self._preinitialize()
     fields = self._initialize(mesh, integrator)
 
-    self.assertAlmostEqual(1.0, integrator.stableTimeStep(mesh)/1.0e+30, places=5)
+    from pylith.utils.utils import maxscalar
+    self.assertAlmostEqual(1.0, integrator.stableTimeStep(mesh)/maxscalar(), 7)
     return
 
 

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -105,7 +105,8 @@
     (mesh, integrator) = self._preinitialize()
     fields = self._initialize(mesh, integrator)
 
-    self.assertAlmostEqual(1.0, integrator.stableTimeStep(mesh)/1.0e+30, places=5)
+    from pylith.utils.utils import maxscalar
+    self.assertAlmostEqual(1.0, integrator.stableTimeStep(mesh)/maxscalar(), 7)
     return
 
 

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicitLgDeform.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicitLgDeform.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicitLgDeform.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -105,7 +105,8 @@
     (mesh, integrator) = self._preinitialize()
     fields = self._initialize(mesh, integrator)
 
-    self.assertAlmostEqual(1.0, integrator.stableTimeStep(mesh)/1.0e+30, places=5)
+    from pylith.utils.utils import maxscalar
+    self.assertAlmostEqual(1.0, integrator.stableTimeStep(mesh)/maxscalar(), 7)
     return
 
 

Modified: short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticIsotropic3D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticIsotropic3D.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticIsotropic3D.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -85,7 +85,8 @@
     from pylith.topology.Mesh import Mesh
     mesh = Mesh()
     dt = self.material.stableTimeStepImplicit(mesh)
-    self.failIf(dt < 1.0e+30)
+    from pylith.utils.utils import maxdouble
+    self.assertAlmostEqual(1.0, dt/maxdouble())
   
 
   def test_factory(self):

Modified: short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticPlaneStrain.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticPlaneStrain.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticPlaneStrain.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -85,7 +85,8 @@
     from pylith.topology.Mesh import Mesh
     mesh = Mesh()
     dt = self.material.stableTimeStepImplicit(mesh)
-    self.failIf(dt < 1.0e+30)
+    from pylith.utils.utils import maxdouble
+    self.assertAlmostEqual(1.0, dt/maxdouble())
   
 
   def test_factory(self):

Modified: short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticPlaneStress.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticPlaneStress.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticPlaneStress.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -85,7 +85,8 @@
     from pylith.topology.Mesh import Mesh
     mesh = Mesh()
     dt = self.material.stableTimeStepImplicit(mesh)
-    self.failIf(dt < 1.0e+30)
+    from pylith.utils.utils import maxdouble
+    self.assertAlmostEqual(1.0, dt/maxdouble())
   
 
   def test_factory(self):

Modified: short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticStrain1D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticStrain1D.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticStrain1D.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -85,7 +85,8 @@
     from pylith.topology.Mesh import Mesh
     mesh = Mesh()
     dt = self.material.stableTimeStepImplicit(mesh)
-    self.failIf(dt < 1.0e+30)
+    from pylith.utils.utils import maxdouble
+    self.assertAlmostEqual(1.0, dt/maxdouble())
   
 
   def test_factory(self):

Modified: short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticStress1D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticStress1D.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/materials/TestElasticStress1D.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -85,7 +85,8 @@
     from pylith.topology.Mesh import Mesh
     mesh = Mesh()
     dt = self.material.stableTimeStepImplicit(mesh)
-    self.failIf(dt < 1.0e+30)
+    from pylith.utils.utils import maxdouble
+    self.assertAlmostEqual(1.0, dt/maxdouble())
   
 
   def test_factory(self):

Modified: short/3D/PyLith/trunk/unittests/pytests/topology/TestRefineUniform.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/topology/TestRefineUniform.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/topology/TestRefineUniform.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -94,7 +94,7 @@
       fault.inventory.faultLabel = faultGroup
       fault._configure()
 
-      nvertices = fault.numVertices(mesh)
+      nvertices = fault.numVerticesNoMesh(mesh)
       firstFaultVertex = 0
       firstLagrangeVertex = nvertices
       firstFaultCell      = 2*nvertices

Modified: short/3D/PyLith/trunk/unittests/pytests/utils/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/utils/Makefile.am	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/utils/Makefile.am	2012-02-24 01:10:59 UTC (rev 19668)
@@ -28,7 +28,8 @@
 
 noinst_PYTHON = \
 	TestEventLogger.py \
-	TestPetscManager.py
+	TestPetscManager.py \
+	TestConstants.py
 
 
 # End of file 

Copied: short/3D/PyLith/trunk/unittests/pytests/utils/TestConstants.py (from rev 19667, short/3D/PyLith/branches/v1.6-stable/unittests/pytests/utils/TestConstants.py)
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/utils/TestConstants.py	                        (rev 0)
+++ short/3D/PyLith/trunk/unittests/pytests/utils/TestConstants.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2011 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/utils/TestEventLogger.py
+
+## @brief Unit testing of EventLogger object.
+
+import unittest
+
+
+# ----------------------------------------------------------------------
+class TestConstants(unittest.TestCase):
+  """
+  Unit testing of constants.
+  """
+  
+
+  def test_maxdouble(self):
+    """
+    Test maxdouble()
+    """
+    from pylith.utils.utils import maxdouble
+    self.assertAlmostEqual(1.0, maxdouble()/1.0e+99, 7)
+    return
+
+
+  def test_maxflat(self):
+    """
+    Test maxflat()
+    """
+    from pylith.utils.utils import maxfloat
+    self.assertAlmostEqual(1.0, maxfloat()/1.0e+30, 7)
+    return
+
+
+# End of file 

Modified: short/3D/PyLith/trunk/unittests/pytests/utils/testutils.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/utils/testutils.py	2012-02-23 21:59:53 UTC (rev 19667)
+++ short/3D/PyLith/trunk/unittests/pytests/utils/testutils.py	2012-02-24 01:10:59 UTC (rev 19668)
@@ -68,6 +68,9 @@
     from TestEventLogger import TestEventLogger
     suite.addTest(unittest.makeSuite(TestEventLogger))
 
+    from TestConstants import TestConstants
+    suite.addTest(unittest.makeSuite(TestConstants))
+
     return suite
 
 

