[cig-commits] r20815 - short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults

[brad at geodynamics.org]

Author: brad
Date: 2012-10-09 16:53:54 -0700 (Tue, 09 Oct 2012)
New Revision: 20815

Modified:
   short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
Log:
Fixed bug in flipping fault orientation (need to sync decision across processors).

Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2012-10-09 17:41:02 UTC (rev 20814)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2012-10-09 23:53:54 UTC (rev 20815)
@@ -1437,21 +1437,28 @@
     // opposite of what we would want, but we cannot flip the fault
     // normal direction because it is tied to how the cohesive cells
     // are created.
-    assert(vertices->size() > 0);
-    orientationSection->restrictPoint(*vertices->begin(),
-      &orientationVertex[0], orientationVertex.size());
+    int flipLocal = 0;
+    if (vertices->size() > 0) {
+      orientationSection->restrictPoint(*vertices->begin(), &orientationVertex[0], orientationVertex.size());
 
-    assert(2 == spaceDim);
-    const PylithScalar* shearDirVertex = &orientationVertex[0];
-    const PylithScalar* normalDirVertex = &orientationVertex[2];
-    const PylithScalar shearDirDot = 
-      upDir[0] * shearDirVertex[0] + upDir[1] * shearDirVertex[1];
-    const PylithScalar normalDirDot = 
-      upDir[0] * normalDirVertex[0] + upDir[1] * normalDirVertex[1];
+      assert(2 == spaceDim);
+      const PylithScalar* shearDirVertex = &orientationVertex[0];
+      const PylithScalar* normalDirVertex = &orientationVertex[2];
+      const PylithScalar shearDirDot = 
+	upDir[0] * shearDirVertex[0] + upDir[1] * shearDirVertex[1];
+      const PylithScalar normalDirDot = 
+	upDir[0] * normalDirVertex[0] + upDir[1] * normalDirVertex[1];
+      if (normalDirDot * shearDirDot < 0.0) {
+	flipLocal = 1;
+      } // if
+    } // if
+    // Collect flip decisions across all processors
+    int flipGlobal = 0;
+    MPI_Allreduce(&flipLocal, &flipGlobal, 1, MPI_INT, MPI_SUM, faultSieveMesh->comm());
 
     const int ishear = 0;
     const int inormal = 2;
-    if (normalDirDot * shearDirDot < 0.0) {
+    if (flipGlobal > 0) { // flip in any processor wants to flip
       // Flip shear direction
       for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin; 
 	   v_iter != verticesEnd;
@@ -1468,11 +1475,7 @@
       PetscLogFlops(3 + count * 2);
     } // if
 
-  } else if (2 == cohesiveDim && vertices->size() > 0) {
-    // Check orientation of first vertex, if dot product of fault
-    // normal with preferred normal is negative, flip up/down dip
-    // direction.
-    //
+  } else if (2 == cohesiveDim) {
     // 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
@@ -1486,32 +1489,38 @@
     // are used are the opposite of what we would want, but we cannot
     // flip the fault normal direction because it is tied to how the
     // cohesive cells are created.
+    int flipLocal = 0;
+    if (vertices->size() > 0) {
+      orientationSection->restrictPoint(*vertices->begin(), &orientationVertex[0], orientationVertex.size());
 
-    assert(vertices->size() > 0);
-    orientationSection->restrictPoint(*vertices->begin(),
-      &orientationVertex[0], orientationVertex.size());
+      assert(3 == spaceDim);
+      const PylithScalar* dipDirVertex = &orientationVertex[3];
+      const PylithScalar* normalDirVertex = &orientationVertex[6];
+      const PylithScalar dipDirDot = 
+	upDir[0]*dipDirVertex[0] + 
+	upDir[1]*dipDirVertex[1] + 
+	upDir[2]*dipDirVertex[2];
+      const PylithScalar normalDirDot = 
+	upDir[0]*normalDirVertex[0] +
+	upDir[1]*normalDirVertex[1] +
+	upDir[2]*normalDirVertex[2];
 
-    assert(3 == spaceDim);
-    const PylithScalar* dipDirVertex = &orientationVertex[3];
-    const PylithScalar* normalDirVertex = &orientationVertex[6];
-    const PylithScalar dipDirDot = 
-      upDir[0]*dipDirVertex[0] + 
-      upDir[1]*dipDirVertex[1] + 
-      upDir[2]*dipDirVertex[2];
-    const PylithScalar normalDirDot = 
-      upDir[0]*normalDirVertex[0] +
-      upDir[1]*normalDirVertex[1] +
-      upDir[2]*normalDirVertex[2];
+      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).
+	flipLocal = 1;
+      } // if
+    } // if
 
+    // Collect flip decisions across all processors
+    int flipGlobal = 0;
+    MPI_Allreduce(&flipLocal, &flipGlobal, 1, MPI_INT, MPI_SUM, faultSieveMesh->comm());
+
     const int istrike = 0;
     const int idip = 3;
     const int inormal = 6;
-    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).
-
+    if (flipGlobal > 0) { // flip in any processor wants to flip
       // Flip dip direction
       for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
 	     != verticesEnd; ++v_iter) {