[cig-commits] r18795 - short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults

[brad at geodynamics.org]

Author: brad
Date: 2011-07-22 08:22:56 -0700 (Fri, 22 Jul 2011)
New Revision: 18795

Modified:
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
Log:
Fixed bug in friction implementation. Account for possibility that shear traction may be zero (opening + sliding).

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc	2011-07-22 01:55:47 UTC (rev 18794)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc	2011-07-22 15:22:56 UTC (rev 18795)
@@ -493,6 +493,8 @@
     dLagrangeTpdtSection->updatePoint(v_fault, &dLagrangeTpdtVertex[0]);
   } // for
 
+  dLagrangeTpdtSection->view("dLagrangeTpdt");
+
   // Solve sensitivity problem for negative side of the fault.
   bool negativeSide = true;
   _sensitivityUpdateJacobian(negativeSide, jacobian, *fields);
@@ -1965,13 +1967,18 @@
     if (tractionShearMag > frictionStress || slipRateMag > 0.0) {
       // traction is limited by friction, so have sliding OR
       // friction exceeds traction due to overshoot in slip
-      
-      // Update traction increment based on value required to stick
-      // versus friction
-      const double dlp = -(tractionShearMag - frictionStress) * area *
-	tractionTpdt[0] / tractionShearMag;
-      (*dLagrangeTpdt)[0] = dlp;
-      (*dLagrangeTpdt)[1] = 0.0;
+
+      if (tractionShearMag > 0.0) {
+	// Update traction increment based on value required to stick
+	// versus friction
+	const double dlp = -(tractionShearMag - frictionStress) * area *
+	  tractionTpdt[0] / tractionShearMag;
+	(*dLagrangeTpdt)[0] = dlp;
+	(*dLagrangeTpdt)[1] = 0.0;
+      } else {
+	(*dLagrangeTpdt)[0] = -(*dLagrangeTpdt)[0];
+	(*dLagrangeTpdt)[1] = 0.0;
+      } // if/else
     } else {
       // friction exceeds value necessary to stick
       // no changes to solution
@@ -2015,16 +2022,22 @@
       // traction is limited by friction, so have sliding OR
       // friction exceeds traction due to overshoot in slip
       
-      // Update traction increment based on value required to stick
-      // versus friction
-      const double dlp = -(tractionShearMag - frictionStress) * area *
-	tractionTpdt[0] / tractionShearMag;
-      const double dlq = -(tractionShearMag - frictionStress) * area *
-	tractionTpdt[1] / tractionShearMag;
-
-      (*dLagrangeTpdt)[0] = dlp;
-      (*dLagrangeTpdt)[1] = dlq;
-      (*dLagrangeTpdt)[2] = 0.0;
+      if (tractionShearMag > 0.0) {
+	// Update traction increment based on value required to stick
+	// versus friction
+	const double dlp = -(tractionShearMag - frictionStress) * area *
+	  tractionTpdt[0] / tractionShearMag;
+	const double dlq = -(tractionShearMag - frictionStress) * area *
+	  tractionTpdt[1] / tractionShearMag;
+	
+	(*dLagrangeTpdt)[0] = dlp;
+	(*dLagrangeTpdt)[1] = dlq;
+	(*dLagrangeTpdt)[2] = 0.0;
+      } else {
+	(*dLagrangeTpdt)[0] = -(*dLagrangeTpdt)[0];
+	(*dLagrangeTpdt)[0] = -(*dLagrangeTpdt)[0];
+	(*dLagrangeTpdt)[2] = 0.0;
+      } // if/else	
       
     } else {
       // else friction exceeds value necessary, so stick