[cig-commits] r16972 - short/3D/PyLith/trunk/libsrc/faults

[brad at geodynamics.org]

Author: brad
Date: 2010-06-09 22:17:41 -0700 (Wed, 09 Jun 2010)
New Revision: 16972

Modified:
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
Log:
Removed ambiguity in dip-slip for horizontal fault.

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc	2010-06-10 04:42:34 UTC (rev 16971)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc	2010-06-10 05:17:41 UTC (rev 16972)
@@ -1717,10 +1717,12 @@
     // normal with preferred normal is negative, flip up/down dip
     // direction.
     //
-    // If the user gives the correct normal direction (points from
-    // footwall to ahanging wall), we should end up with
-    // left-lateral-slip, reverse-slip, and fault-opening for positive
-    // slip values.
+    // Check orientation of first vertex, (1) if dot product of the
+    // normal-dir with preferred up-dir is positive, then we want dot
+    // product of shear-dir and preferred up-dir to be positive and
+    // (2) if the dot product of the normal-dir with preferred up-dir
+    // is negative, then we want the dot product of the shear-dir and
+    // preferred up-dir to be negative.
     //
     // When we flip the up/down dip direction, we create a left-handed
     // strike/dip/normal coordinate system, but it gives the correct
@@ -1748,16 +1750,21 @@
     const int istrike = 0;
     const int idip = 3;
     const int inormal = 6;
-    if (dipDirDot * normalDirDot < 0.0) {
+    if (dipDirDot * normalDirDot < 0.0 ||
+	fabs(normalDirVertex[2] + 1.0) < 0.001) {
+      // if fault normal is (0,0,+-1) then up-dir dictates reverse
+      // motion for case with normal (0,0,1), so we reverse the dip-dir
+      // if we have (0,0,-1).
+
       // Flip dip direction
       for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
-          != verticesEnd; ++v_iter) {
+	     != verticesEnd; ++v_iter) {
         orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
-          orientationVertex.size());
+					  orientationVertex.size());
         assert(9 == orientationSection->getFiberDimension(*v_iter));
         for (int iDim = 0; iDim < 3; ++iDim) // flip dip
           orientationVertex[idip + iDim] *= -1.0;
-
+	
         // Update direction
         orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
       } // for