[cig-commits] r15948 - in short/3D/PyLith/branches/pylith-friction/libsrc: faults problems

[brad at geodynamics.org]

Author: brad
Date: 2009-11-10 17:58:10 -0800 (Tue, 10 Nov 2009)
New Revision: 15948

Modified:
   short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc
   short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc
Log:
Fixed bug in friction implementation. Need to update PetscVec of solution after constraining Lagrange multiplier values.

Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc	2009-11-11 00:44:18 UTC (rev 15947)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc	2009-11-11 01:58:10 UTC (rev 15948)
@@ -800,7 +800,8 @@
 	    const double Spq = Cpx*Cqx*Aixjx + Cpy*Cqy*Aiyjy;
 	    const double Sqq = Cqx*Cqx*Aixjx + Cqy*Cqy*Aiyjy;
 
-	    if (tractionTpdtVertex[1]+tractionInitialVertex[1] < 0) {
+	    if (tractionTpdtVertex[1]+tractionInitialVertex[1] < 0 &&
+		slipVertex[1] == 0) {
 	      // if in compression
 	      std::cout << "FAULT IN COMPRESSION" << std::endl;
 	      const double friction =

Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc	2009-11-11 00:44:18 UTC (rev 15947)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc	2009-11-11 01:58:10 UTC (rev 15948)
@@ -153,10 +153,14 @@
     _submeshIntegrators[i]->timeStep(_dt);
     _submeshIntegrators[i]->constrainSolnSpace(_fields, _t, *_jacobian);
   } // for
+
+  if (0 != tmpSolutionVec) {
+    topology::Field<topology::Mesh>& solution = _fields->solution();
+    solution.scatterSectionToVector(*tmpSolutionVec);
+  } // if
   // :KLUDGE: END
   // ===================================================================
 
-
   // Set residual to zero.
   topology::Field<topology::Mesh>& residual = _fields->get("residual");
   residual.zero();