[cig-commits] [commit] knepley/fix-faults-parallel: Fault: We need a consistent closure for clamped cohesive cells - We get the two end faces of the prism and concatenate their closure (0f8af6e)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/pylith

On branch  : knepley/fix-faults-parallel
Link       : https://github.com/geodynamics/pylith/compare/eeb238185ffa74426d985c5e35091fd1d3ad8db4...0f8af6e7a83ed15f72a991e2a2f8b66fd29f9c12

>---------------------------------------------------------------

commit 0f8af6e7a83ed15f72a991e2a2f8b66fd29f9c12
Author: Matthew G. Knepley <knepley at gmail.com>
Date:   Tue May 27 18:54:12 2014 -0500

    Fault: We need a consistent closure for clamped cohesive cells
    - We get the two end faces of the prism and concatenate their closure


>---------------------------------------------------------------

0f8af6e7a83ed15f72a991e2a2f8b66fd29f9c12
 libsrc/pylith/faults/FaultCohesiveDyn.cc | 40 +++++++++++++++++++++++---------
 1 file changed, 29 insertions(+), 11 deletions(-)

diff --git a/libsrc/pylith/faults/FaultCohesiveDyn.cc b/libsrc/pylith/faults/FaultCohesiveDyn.cc
index 95be35a..7a2676a 100644
--- a/libsrc/pylith/faults/FaultCohesiveDyn.cc
+++ b/libsrc/pylith/faults/FaultCohesiveDyn.cc
@@ -1746,25 +1746,42 @@ pylith::faults::FaultCohesiveDyn::_sensitivityUpdateJacobian(const bool negative
   int_array indicesPerm(subnrows);
   for (PetscInt c = 0; c < numCohesiveCells; ++c) {
     // Get cone for cohesive cell
-    PetscInt* closure = NULL;
-    PetscInt closureSize, q = 0;
-    err = DMPlexGetTransitiveClosure(dmMesh, cellsCohesive[c], PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
-
+    const PetscInt *cone;
+    PetscInt        coneSize;
+    PetscInt       *closureA = NULL, *closureB = NULL;
+    PetscInt        closureSizeA, closureSizeB, q;
+
+    err = DMPlexGetCone(dmMesh, cellsCohesive[c], &cone);PYLITH_CHECK_ERROR(err);
+    err = DMPlexGetConeSize(dmMesh, cellsCohesive[c], &coneSize);PYLITH_CHECK_ERROR(err);
+    assert(coneSize >= 4);
+    err = DMPlexGetTransitiveClosure(dmMesh, cone[0], PETSC_TRUE, &closureSizeA, &closureA);PYLITH_CHECK_ERROR(err);
+    // Filter out non-vertices
+    q = 0;
+    for(PetscInt p = 0; p < closureSizeA*2; p += 2) {
+      if ((closureA[p] >= vStart) && (closureA[p] < vEnd)) {
+        closureA[q] = closureA[p];
+        ++q;
+      } // if
+    } // for
+    closureSizeA = q;
+    err = DMPlexGetTransitiveClosure(dmMesh, cone[1], PETSC_TRUE, &closureSizeB, &closureB);PYLITH_CHECK_ERROR(err);
     // Filter out non-vertices
-    for (PetscInt p = 0; p < closureSize*2; p += 2) {
-      if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
-        closure[q] = closure[p];
+    q = 0;
+    for(PetscInt p = 0; p < closureSizeB*2; p += 2) {
+      if ((closureB[p] >= vStart) && (closureB[p] < vEnd)) {
+        closureB[q] = closureB[p];
         ++q;
       } // if
     } // for
-    closureSize = q;
-    assert(closureSize == 2*numBasis);
+    closureSizeB = q;
+    assert(closureSizeA == numBasis);
+    assert(closureSizeB == numBasis);
 
     // Get indices
     for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
       // negative side of the fault: iCone=0
       // positive side of the fault: iCone=1
-      const int v_domain = closure[iCone*numBasis+iBasis];
+      const int v_domain = iCone ? closureB[iBasis] :  closureA[iBasis];
       PetscInt goff;
 
       err = PetscSectionGetOffset(solutionDomainGlobalSection, v_domain, &goff);PYLITH_CHECK_ERROR(err);
@@ -1773,7 +1790,8 @@ pylith::faults::FaultCohesiveDyn::_sensitivityUpdateJacobian(const bool negative
       } // for
 
     } // for
-    err = DMPlexRestoreTransitiveClosure(dmMesh, cellsCohesive[c], PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
+    err = DMPlexRestoreTransitiveClosure(dmMesh, cone[0], PETSC_TRUE, &closureSizeA, &closureA);PYLITH_CHECK_ERROR(err);
+    err = DMPlexRestoreTransitiveClosure(dmMesh, cone[1], PETSC_TRUE, &closureSizeB, &closureB);PYLITH_CHECK_ERROR(err);
 
     for (int i=0; i < subnrows; ++i) {
       indicesPerm[i]  = i;