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

[brad at geodynamics.org]

Author: brad
Date: 2011-08-19 09:47:18 -0700 (Fri, 19 Aug 2011)
New Revision: 18838

Modified:
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
Log:
Fixed bug in adjusting slip in quasistatic solve to handle changes in fault opening displacement. Must account for current slip in addition to slip increment.

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-08-19 16:43:36 UTC (rev 18837)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc	2011-08-19 16:47:18 UTC (rev 18838)
@@ -521,6 +521,9 @@
     orientationSection->restrictPoint(v_fault, &orientationVertex[0],
     orientationVertex.size());
 
+    // Get slip
+    slipSection->restrictPoint(v_fault, &slipVertex[0], slipVertex.size());
+
     // Get change in relative displacement.
     const double* dispRelVertex = dispRelSection->restrictPoint(v_fault);
     assert(0 != dispRelVertex);
@@ -555,12 +558,12 @@
     // Do not allow fault interpenetration and set fault opening to
     // zero if fault is under compression.
     const int indexN = spaceDim - 1;
-    if (dSlipVertex[indexN] < 0.0)
-      dSlipVertex[indexN] = 0.0;
+    if (slipVertex[indexN] + dSlipVertex[indexN] < 0.0)
+      dSlipVertex[indexN] = -dSlipVertex[indexN];
     const double lagrangeTpdtNormal = lagrangeTVertex[indexN] + 
       lagrangeTIncrVertex[indexN] + dLagrangeTpdtVertex[indexN];
     if (lagrangeTpdtNormal < 0.0)
-      dSlipVertex[indexN] = 0.0;
+      dSlipVertex[indexN] = -slipVertex[indexN];
 
     // Set change in slip.
     assert(dSlipVertex.size() ==
@@ -1792,7 +1795,10 @@
   // Update section view of field.
   solution.scatterVectorToSection();
 
-  //solution.view("SENSITIVITY SOLUTION"); // DEBUGGING
+#if 0 // DEBUGGING
+  residual.view("SENSITIVITY RESIDUAL");
+  solution.view("SENSITIVITY SOLUTION");
+#endif
 } // _sensitivitySolve
 
 // ----------------------------------------------------------------------