[cig-commits] [commit] knepley/upgrade-petsc-interface: Updated fault C++ unit tests. (dcfe4fb)

cig_noreply at geodynamics.org cig_noreply at geodynamics.org
Fri Nov 8 17:07:20 PST 2013


Repository : ssh://geoshell/pylith

On branch  : knepley/upgrade-petsc-interface
Link       : https://github.com/geodynamics/pylith/compare/44a899ad4662f73e26966c1cf40ca2eb5790ee6e...1cae3da52b30afc7095c28b401d4e90b5f51a0f4

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

commit dcfe4fbef7f58be6a7a12af55e76a1f7b29f8c04
Author: Brad Aagaard <baagaard at usgs.gov>
Date:   Fri Nov 8 16:48:49 2013 -0800

    Updated fault C++ unit tests.
    
    Updated for interpolated meshes. Some fixes due to corruption of test
    data in previous update for interpolated meshes.


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

dcfe4fbef7f58be6a7a12af55e76a1f7b29f8c04
 libsrc/pylith/faults/FaultCohesiveDyn.cc           |  33 +--
 libsrc/pylith/faults/FaultCohesiveLagrange.cc      |   2 +-
 unittests/libtests/faults/TestFaultCohesiveDyn.cc  |   6 +-
 .../libtests/faults/TestFaultCohesiveImpulses.cc   |   4 +-
 .../libtests/faults/data/CohesiveDynDataHex8.cc    |  48 ++--
 .../libtests/faults/data/CohesiveDynDataQuad4.cc   |  24 +-
 .../libtests/faults/data/CohesiveDynDataTet4.cc    |  66 +++---
 .../libtests/faults/data/CohesiveDynDataTri3.cc    |  22 +-
 .../libtests/faults/data/CohesiveDynDataTri3d.cc   |  60 ++---
 unittests/libtests/faults/data/cohesivedyn.py      | 258 +++++++++++----------
 10 files changed, 249 insertions(+), 274 deletions(-)

diff --git a/libsrc/pylith/faults/FaultCohesiveDyn.cc b/libsrc/pylith/faults/FaultCohesiveDyn.cc
index 1be2dc8..3afb10a 100644
--- a/libsrc/pylith/faults/FaultCohesiveDyn.cc
+++ b/libsrc/pylith/faults/FaultCohesiveDyn.cc
@@ -963,7 +963,7 @@ pylith::faults::FaultCohesiveDyn::constrainSolnSpace(topology::SolutionFields* c
       dDispTIncrVertexP[iDim] = +0.5*dDispRelVertex[iDim];
     } // for
 
-#if 1 // debugging
+#if 0 // debugging
     std::cout << "v_fault: " << v_fault;
     std::cout << ", tractionTpdtVertex: ";
     for (int iDim=0; iDim < spaceDim; ++iDim)
@@ -1747,7 +1747,6 @@ pylith::faults::FaultCohesiveDyn::_sensitivityUpdateJacobian(const bool negative
         indicesGlobal[iB+iDim] = gind + iDim;
       } // for
 
-      std::cout << "v_domain["<<iBasis<<"]: " << v_domain << ", globalIndex: " << goff << std::endl;
     } // for
     err = DMPlexRestoreTransitiveClosure(dmMesh, cellsCohesive[c], PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
 
@@ -1762,45 +1761,19 @@ pylith::faults::FaultCohesiveDyn::_sensitivityUpdateJacobian(const bool negative
     cellsIS[c] = NULL;
     err = ISCreateGeneral(PETSC_COMM_SELF, indicesGlobal.size(), &indicesGlobal[0], PETSC_COPY_VALUES, &cellsIS[c]);PYLITH_CHECK_ERROR(err);
 
-#if 1 // DEBUGGING
-    std::cout << "indicesGlobal:";
-    for (int ii=0; ii < indicesGlobal.size(); ++ii) {
-      std::cout << " " << indicesGlobal[ii];
-    }
-    std::cout << std::endl;
-
-    std::cout << "indicesLocal:";
-    for (int ii=0; ii < subnrows; ++ii) {
-      std::cout << " " << indicesLocal[c*subnrows+ii];
-    }
-    std::cout << std::endl;
-
-#endif
-
   } // for
 
   PetscMat* submatrices = NULL;
   err = MatGetSubMatrices(jacobianDomainMatrix, numCohesiveCells, cellsIS, cellsIS, MAT_INITIAL_MATRIX, &submatrices);PYLITH_CHECK_ERROR(err);
 
-  _faultMesh->view("fault", "::ascii_info_detail");
-
   for (PetscInt c = 0; c < numCohesiveCells; ++c) {
     // Get values for submatrix associated with cohesive cell
     jacobianSubCell = 0.0;
     err = MatGetValues(submatrices[c], subnrows, &indicesLocal[c*subnrows], subnrows, &indicesLocal[c*subnrows],
                        &jacobianSubCell[0]);PYLITH_CHECK_ERROR_MSG(err, "Restrict from PETSc Mat failed.");
 
-    std::cout << "SUBMATRIX " << c << std::endl;
-    MatView(submatrices[c], PETSC_VIEWER_STDOUT_WORLD);
-    std::cout << "jacobianSubCell: ";
-    for (int ii=0; ii < subnrows*subnrows; ++ii) {
-      std::cout << " " << jacobianSubCell[ii];
-    }
-    std::cout << std::endl;
-
     // Insert cell contribution into PETSc Matrix
     PetscInt c_fault = _cohesiveToFault[cellsCohesive[c]];
-    std::cout << "c_fault: " << c_fault << std::endl;
 
     err = DMPlexMatSetClosure(faultDMMesh, solutionFaultSection, solutionFaultGlobalSection,  jacobianFaultMatrix, c_fault, &jacobianSubCell[0], INSERT_VALUES);PYLITH_CHECK_ERROR_MSG(err, "Update to PETSc Mat failed.");
 
@@ -1813,7 +1786,7 @@ pylith::faults::FaultCohesiveDyn::_sensitivityUpdateJacobian(const bool negative
 
   _jacobian->assemble("final_assembly");
 
-#if 1 // DEBUGGING
+#if 0 // DEBUGGING
   //std::cout << "DOMAIN JACOBIAN" << std::endl;
   //jacobian.view();
   std::cout << "SENSITIVITY JACOBIAN" << std::endl;
@@ -1944,7 +1917,7 @@ pylith::faults::FaultCohesiveDyn::_sensitivitySolve(void)
   // Update section view of field.
   solution.scatterGlobalToLocal();
 
-#if 1 // DEBUGGING
+#if 0 // DEBUGGING
   residual.view("SENSITIVITY RESIDUAL");
   solution.view("SENSITIVITY SOLUTION");
 #endif
diff --git a/libsrc/pylith/faults/FaultCohesiveLagrange.cc b/libsrc/pylith/faults/FaultCohesiveLagrange.cc
index 1f7417e..cb318f6 100644
--- a/libsrc/pylith/faults/FaultCohesiveLagrange.cc
+++ b/libsrc/pylith/faults/FaultCohesiveLagrange.cc
@@ -1111,7 +1111,7 @@ void pylith::faults::FaultCohesiveLagrange::_initializeCohesiveInfo(const topolo
         _cohesiveVertices[index].positive = v_positive;
         _cohesiveVertices[index].negative = v_negative;
         _cohesiveVertices[index].fault = v_fault;
-#if 1
+#if 0 // DEBUGGING
         std::cout << "cohesiveVertices[" << index << "]: "
 		  << "l: " << e_lagrange
 		  << ", p: " << v_positive
diff --git a/unittests/libtests/faults/TestFaultCohesiveDyn.cc b/unittests/libtests/faults/TestFaultCohesiveDyn.cc
index db97d5a..80f8c1a 100644
--- a/unittests/libtests/faults/TestFaultCohesiveDyn.cc
+++ b/unittests/libtests/faults/TestFaultCohesiveDyn.cc
@@ -178,7 +178,7 @@ pylith::faults::TestFaultCohesiveDyn::testInitialize(void)
   CPPUNIT_ASSERT_EQUAL(_data->numConstraintEdges, vEnd-vStart);
 
   // Check orientation
-  fault._fields->get("orientation").view("ORIENTATION"); // DEBUGGING
+  //fault._fields->get("orientation").view("ORIENTATION"); // DEBUGGING
   topology::VecVisitorMesh orientationVisitor(fault._fields->get("orientation"));
   const PetscScalar* orientationArray = orientationVisitor.localArray();
 
@@ -333,7 +333,7 @@ pylith::faults::TestFaultCohesiveDyn::testConstrainSolnSpaceSlip(void)
 
   fault.updateStateVars(t, &fields);
 
-  solution.view("SOLUTION"); // DEBUGGING
+  //solution.view("SOLUTION"); // DEBUGGING
 
   { // Check solution values
     // Lagrange multipliers should be adjusted according to friction
@@ -425,7 +425,7 @@ pylith::faults::TestFaultCohesiveDyn::testConstrainSolnSpaceOpen(void)
 
   fault.updateStateVars(t, &fields);
 
-  solution.view("SOLUTION"); // DEBUGGING
+  //solution.view("SOLUTION"); // DEBUGGING
 
   { // Check solution values
     // Lagrange multipliers should be set to zero as reflected in the
diff --git a/unittests/libtests/faults/TestFaultCohesiveImpulses.cc b/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
index 026e983..1522bce 100644
--- a/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
+++ b/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
@@ -196,7 +196,7 @@ pylith::faults::TestFaultCohesiveImpulses::testInitialize(void)
   CPPUNIT_ASSERT_EQUAL(_data->numConstraintEdges, vEnd-vStart);
 
   // Check orientation
-  fault._fields->get("orientation").view("ORIENTATION"); // DEBUGGING
+  //fault._fields->get("orientation").view("ORIENTATION"); // DEBUGGING
 
   topology::VecVisitorMesh orientationVisitor(fault._fields->get("orientation"));
   const PetscScalar* orientationArray = orientationVisitor.localArray();
@@ -279,7 +279,7 @@ pylith::faults::TestFaultCohesiveImpulses::testIntegrateResidual(void)
   topology::Field& residual = fields.get("residual");
   residual.zero();
   fault.integrateResidual(residual, t, &fields);
-  residual.view("RESIDUAL"); // DEBUGGING
+  //residual.view("RESIDUAL"); // DEBUGGING
 
   topology::VecVisitorMesh residualVisitor(residual);
   const PetscScalar* residualArray = residualVisitor.localArray();CPPUNIT_ASSERT(residualArray);
diff --git a/unittests/libtests/faults/data/CohesiveDynDataHex8.cc b/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
index 6fe5995..d566093 100644
--- a/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
+++ b/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
@@ -1423,18 +1423,18 @@ const PylithScalar pylith::faults::CohesiveDynDataHex8::_fieldIncrSlipE[20*3] =
    1.200000000000,   2.200000000000,   0.200000000000,
    1.300000000000,   2.300000000000,   0.300000000000,
    1.400000000000,   2.400000000000,   0.400000000000,
-   1.500000000000,   2.500001237431,   0.500000808427,
-   1.600000000000,   2.600000916263,   0.600000905670,
-   1.700000000000,   2.700000658910,   0.700000753748,
-   1.800000000000,   2.800000898799,   0.800000868468,
+   1.500000000000,   2.500001245294,   0.500000812227,
+   1.600000000000,   2.600000694980,   0.600000786664,
+   1.700000000000,   2.700000891763,   0.700000882544,
+   1.800000000000,   2.800000903108,   0.800000871357,
    1.900000000000,   2.900000000000,   0.900000000000,
    1.000000000000,   2.000000000000,   0.000000000000,
    1.100000000000,   2.100000000000,   0.100000000000,
    1.200000000000,   2.200000000000,   0.200000000000,
-   1.500000000000,   2.499998762569,   0.499999191573,
-   1.600000000000,   2.599999083737,   0.599999094330,
-   1.700000000000,   2.699999341090,   0.699999246252,
-   1.800000000000,   2.799999101201,   0.799999131532,
+   1.500000000000,   2.499998754706,   0.499999187773,
+   1.600000000000,   2.599999305020,   0.599999213336,
+   1.700000000000,   2.699999108237,   0.699999117456,
+   1.800000000000,   2.799999096892,   0.799999128643,
   -1.400000000000,   0.328479469127,  -1.239953753608,
   -1.600000000000,   0.293941190358,  -1.262585961634,
   -1.800000000000,   0.259995967920,  -1.286431883494,
@@ -1442,10 +1442,10 @@ const PylithScalar pylith::faults::CohesiveDynDataHex8::_fieldIncrSlipE[20*3] =
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataHex8::_slipSlipE[] = {
-   0.699997525137,   0.799998383146,   0.000000000000,
-   0.899998167474,   0.999998188660,   0.000000000000,
-   0.999998682181,   0.899998492503,   0.000000000000,
-   0.799998202402,   0.699998263064,   0.000000000000,
+   0.699997509413,   0.799998375546,   0.000000000000,
+   0.899998610039,   0.999998426672,   0.000000000000,
+   0.999998216473,   0.899998234913,   0.000000000000,
+   0.799998193784,   0.699998257287,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -1481,18 +1481,18 @@ const PylithScalar pylith::faults::CohesiveDynDataHex8::_fieldIncrOpenE[] = {
    1.200000000000,   2.200000000000,   0.200000000000,
    1.300000000000,   2.300000000000,   0.300000000000,
    1.400000000000,   2.400000000000,   0.400000000000,
-   1.500000000000,   2.500007868394,   0.500005286699,
-   1.600000000000,   2.600006129978,   0.600006000739,
-   1.700000000000,   2.700004813062,   0.700005231911,
-   1.800000000000,   2.800006278922,   0.800006060966,
+   1.500000000000,   2.500007882076,   0.500005293315,
+   1.600000000000,   2.600004874790,   0.600005290605,
+   1.700000000000,   2.700006088248,   0.700005959177,
+   1.800000000000,   2.800006286383,   0.800006065901,
    1.900000000000,   2.900000000000,   0.900000000000,
    1.000000000000,   2.000000000000,   0.000000000000,
    1.100000000000,   2.100000000000,   0.100000000000,
    1.200000000000,   2.200000000000,   0.200000000000,
-   1.500000000000,   2.499992131606,   0.499994713301,
-   1.600000000000,   2.599993870022,   0.599993999261,
-   1.700000000000,   2.699995186938,   0.699994768089,
-   1.800000000000,   2.799993721078,   0.799993939034,
+   1.500000000000,   2.499992117924,   0.499994706685,
+   1.600000000000,   2.599995125210,   0.599994709395,
+   1.700000000000,   2.699993911752,   0.699994040823,
+   1.800000000000,   2.799993713617,   0.799993934099,
    4.400000000000,  -2.400000000000,  -3.400000000000,
    4.600000000000,  -2.600000000000,  -3.600000000000,
    4.800000000000,  -2.800000000000,  -3.800000000000,
@@ -1500,10 +1500,10 @@ const PylithScalar pylith::faults::CohesiveDynDataHex8::_fieldIncrOpenE[] = {
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataHex8::_slipOpenE[] = {
-   0.699984263212,   0.799989426603,   0.000000000000,
-   0.899987740045,   0.999987998523,   0.000000000000,
-   0.999990373875,   0.899989536177,   0.000000000000,
-   0.799987442156,   0.699987878069,   0.000000000000,
+   0.699984235848,   0.799989413370,   0.000000000000,
+   0.899990250419,   0.999989418790,   0.000000000000,
+   0.999987823504,   0.899988081646,   0.000000000000,
+   0.799987427233,   0.699987868197,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
diff --git a/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc b/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
index ab9c989..600b842 100644
--- a/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
+++ b/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
@@ -339,19 +339,19 @@ const PylithScalar pylith::faults::CohesiveDynDataQuad4::_fieldIncrSlip[] = {
 const PylithScalar pylith::faults::CohesiveDynDataQuad4::_fieldIncrSlipE[] = {
    1.100000000000,   2.100000000000,
    1.200000000000,   2.200000000000,
-   1.200000000000,   2.200960513755,
-   1.300000000000,   2.300184361537,
+   1.200000000000,   2.200157498929,
+   1.300000000000,   2.300968885890,
    1.500000000000,   2.500000000000,
    1.600000000000,   2.600000000000,
-   1.200000000000,   2.199039486245,
-   1.300000000000,   2.299815638463,
+   1.200000000000,   2.199842501071,
+   1.300000000000,   2.299031114110,
   -1.600000000000,  -3.480000000000,
   -1.800000000000,  -3.440000000000,
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataQuad4::_slipSlipE[] = {
-   0.398078972490,   0.000000000000,
-   0.499631276926,   0.000000000000,
+   0.399685002143,   0.000000000000,
+   0.498062228219,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -375,19 +375,19 @@ const PylithScalar pylith::faults::CohesiveDynDataQuad4::_fieldIncrOpen[] = {
 const PylithScalar pylith::faults::CohesiveDynDataQuad4::_fieldIncrOpenE[] = {
    1.100000000000,   2.100000000000,
    1.200000000000,   2.200000000000,
-   1.199417086548,   2.200039516147,
-   1.299388886535,   2.300357375764,
+   1.199388886535,   2.200357375764,
+   1.299417086548,   2.300039516147,
    1.500000000000,   2.500000000000,
    1.600000000000,   2.600000000000,
-   1.200582913452,   2.199960483853,
-   1.300611113465,   2.299642624236,
+   1.200611113465,   2.199642624236,
+   1.300582913452,   2.299960483853,
    8.600000000000,  -9.600000000000,
    8.800000000000,  -9.800000000000,
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataQuad4::_slipOpenE[] = {
-   0.399920967705,   0.001165826903,
-   0.499285248472,   0.001222226930,
+   0.399285248472,   0.001222226930,
+   0.499920967705,   0.001165826903,
 };
 
 // ----------------------------------------------------------------------
diff --git a/unittests/libtests/faults/data/CohesiveDynDataTet4.cc b/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
index d0b872e..95d3edd 100644
--- a/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
+++ b/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
@@ -97,9 +97,9 @@ const PylithScalar pylith::faults::CohesiveDynDataTet4::_fieldT[] = {
   7.2, 8.2, 9.2, // 8
   7.3, 8.3, 9.3, // 9
   7.4, 8.4, 9.4, // 10
- -7.7, 8.7, 9.7, // 34
- -7.9, 8.9, 9.9, // 35
- -7.1, 8.1, 9.1, // 36
+ -7.7, 18.7, 19.7, // 34
+ -7.9, 18.9, 19.9, // 35
+ -7.1, 18.1, 19.1, // 36
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataTet4::_jacobian[] = {
@@ -513,9 +513,9 @@ const PylithScalar pylith::faults::CohesiveDynDataTet4::_fieldIncrStick[] = {
   1.2, 2.2, 3.2, // 8
   1.3, 2.3, 3.3, // 9
   1.4, 2.4, 3.4, // 10
- -41.7, 2.7, 3.7, // 34
- -41.9, 2.9, 3.9, // 35
- -41.1, 2.1, 3.1, // 36
+ -81.7, 2.7, 3.7, // 34
+ -81.9, 2.9, 3.9, // 35
+ -81.1, 2.1, 3.1, // 36
 };
 
 // No slip
@@ -538,30 +538,30 @@ const PylithScalar pylith::faults::CohesiveDynDataTet4::_fieldIncrSlip[] = {
   1.2, 2.2, 3.2, // 7
   1.3, 2.3, 3.3, // 8
   1.4, 2.4, 3.4, // 10
- -4.7, 9.7, 3.7, // 34
- -4.9, 9.9, 3.9, // 35
- -4.1, 9.1, 3.1, // 36
+ -4.7, 5.7, 6.7, // 34
+ -4.9, 5.9, 6.9, // 35
+ -4.1, 5.1, 6.1, // 36
 };
 
 // Output
 const PylithScalar pylith::faults::CohesiveDynDataTet4::_fieldIncrSlipE[] = {
    1.100000000000,   2.100000000000,   3.100000000000,
-   1.200000000000,   2.200000405670,   3.200000218150,
-   1.300000000000,   2.300000107994,   3.300000286268,
-   1.400000000000,   2.399999892799,   3.400000205817,
+   1.200000000000,   2.200000916680,   3.200000185945,
+   1.300000000000,   2.299999954294,   3.300000084646,
+   1.400000000000,   2.400000465197,   3.400000369651,
    1.500000000000,   2.500000000000,   3.500000000000,
-   1.200000000000,   2.199999594330,   3.199999781850,
-   1.300000000000,   2.299999892006,   3.299999713732,
-   1.400000000000,   2.400000107201,   3.399999794183,
-  -4.700000000000,  -2.685831659394,  -5.320116534558,
-  -4.900000000000,  -2.708904191586,  -5.355472225738,
-  -4.100000000000,  -2.618821803722,  -5.212187558454,
+   1.200000000000,   2.199999083320,   3.199999814055,
+   1.300000000000,   2.300000045706,   3.299999915354,
+   1.400000000000,   2.399999534803,   3.399999630349,
+  -4.700000000000, -13.650158708856, -14.236237291549,
+  -4.900000000000, -13.683824215808, -14.263164878373,
+  -4.100000000000, -13.548480328050, -14.156107942537,
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataTet4::_slipSlipE[] = {
-  -0.000000811341,  -0.000000436300,   0.000000000000,
-  -0.000000215989,  -0.000000572537,   0.000000000000,
-   0.000000214401,  -0.000000411634,   0.000000000000,
+  -0.000001833360,  -0.000000371890,   0.000000000000,
+   0.000000091413,  -0.000000169291,   0.000000000000,
+  -0.000000930394,  -0.000000739303,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -585,22 +585,22 @@ const PylithScalar pylith::faults::CohesiveDynDataTet4::_fieldIncrOpen[] = {
 // Output
 const PylithScalar pylith::faults::CohesiveDynDataTet4::_fieldIncrOpenE[] = {
    1.100000000000,   2.100000000000,   3.100000000000,
-   1.199998457080,   2.200003764775,   3.200002359542,
-   1.299999462980,   2.300002309053,   3.300002198607,
-   1.399997561706,   2.400002145902,   3.400002499006,
+   1.199999161233,   2.200004202086,   3.200002477003,
+   1.299998110594,   2.300001999917,   3.300002482177,
+   1.400000000000,   2.400002588150,   3.400002474071,
    1.500000000000,   2.500000000000,   3.500000000000,
-   1.200001542920,   2.199996235225,   3.199997640458,
-   1.300000537020,   2.299997690947,   3.299997801393,
-   1.400002438294,   2.399997854098,   3.399997500994,
-   7.700000000000,  -8.700000000000,  -9.700000000000,
-   7.900000000000,  -8.900000000000,  -9.900000000000,
-   7.100000000000,  -8.100000000000,  -9.100000000000,
+   1.200000838767,   2.199995797914,   3.199997522997,
+   1.300001889406,   2.299998000083,   3.299997517823,
+   1.400000000000,   2.399997411850,   3.399997525929,
+   7.700000000000, -18.700000000000, -19.700000000000,
+   7.900000000000, -18.900000000000, -19.900000000000,
+   7.100000000000, -18.100000000000, -19.100000000000,
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataTet4::_slipOpenE[] = {
-  -0.000007529551,  -0.000004719084,   0.000003085841,
-  -0.000004618106,  -0.000004397213,   0.000001074039,
-  -0.000004291805,  -0.000004998013,   0.000004876589,
+  -0.000008404171,  -0.000004954006,   0.000001677534,
+  -0.000003999835,  -0.000004964354,   0.000003778811,
+  -0.000005176299,  -0.000004948142,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
diff --git a/unittests/libtests/faults/data/CohesiveDynDataTri3.cc b/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
index 08bcea0..ca2cc76 100644
--- a/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
+++ b/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
@@ -304,18 +304,18 @@ const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrSlip[] = {
 // Output
 const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrSlipE[] = {
    9.100000000000,   7.100000000000,
-   9.200000000000,   7.200993561267,
-   9.300000000000,   7.300159846316,
+   9.200000000000,   7.200190546440,
+   9.300000000000,   7.300983993112,
    9.400000000000,   7.400000000000,
-   9.200000000000,   7.199006438733,
-   9.300000000000,   7.299840153684,
+   9.200000000000,   7.199809453560,
+   9.300000000000,   7.299016006888,
   -1.600000000000,  -3.480000000000,
   -1.800000000000,  -3.440000000000,
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataTri3::_slipSlipE[] = {
-  -0.001987122533,   0.000000000000,
-  -0.000319692633,   0.000000000000,
+  -0.000381092881,   0.000000000000,
+  -0.001967986224,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -336,18 +336,18 @@ const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrOpen[] = {
 // Output
 const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrOpenE[] = {
    9.100000000000,   7.100000000000,
-   9.199943020362,   7.200360547702,
-   9.299826714877,   7.300414257705,
+   9.199826714877,   7.200414257705,
+   9.299943020362,   7.300360547702,
    9.400000000000,   7.400000000000,
-   9.200056979638,   7.199639452298,
-   9.300173285123,   7.299585742295,
+   9.200173285123,   7.199585742295,
+   9.300056979638,   7.299639452298,
    8.600000000000,  -9.600000000000,
    8.800000000000,  -9.800000000000,
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataTri3::_slipOpenE[] = {
-  -0.000721095405,   0.000113959275,
   -0.000828515410,   0.000346570247,
+  -0.000721095405,   0.000113959275,
 };
 
 // ----------------------------------------------------------------------
diff --git a/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc b/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
index 21090ca..3428474 100644
--- a/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
+++ b/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
@@ -487,31 +487,31 @@ const PylithScalar pylith::faults::CohesiveDynDataTri3d::_fieldIncrSlip[] = {
   1.2, 2.2, // 10
   1.3, 2.3, // 11
   1.5, 2.5, // 12
- -1.8,-0.8, // 26
- -1.0, 0.1, // 27
-  1.2,-0.2, // 28
+ -1.8,+3.6, // 26
+ -1.0, 1.1, // 27
+  1.7,-1.2, // 28
 };
 
 // Output
 const PylithScalar pylith::faults::CohesiveDynDataTri3d::_fieldIncrSlipE[] = {
    1.100000000000,   2.100000000000,
-   1.200805761009,   2.199194238991,
-   1.300000000000,   2.299628356338,
+   1.199970441179,   2.200029558821,
+   1.300000000000,   2.299168310929,
    1.400000000000,   2.400000000000,
-   1.499716364403,   2.500000000000,
+   1.499489949674,   2.500000000000,
    1.600000000000,   2.600000000000,
-   1.199194238991,   2.200805761009,
-   1.300000000000,   2.300371643662,
-   1.500283635597,   2.500000000000,
-   4.520000000000,  -1.920000000000,
-   1.000000000000,  -1.600000000000,
-  -0.560000000000,   0.200000000000,
+   1.200029558821,   2.199970441179,
+   1.300000000000,   2.300831689071,
+   1.500510050326,   2.500000000000,
+  -1.640000000000,   3.440000000000,
+  -1.000000000000,  -1.600000000000,
+   0.040000000000,  -1.200000000000,
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataTri3d::_slipSlipE[] = {
-  -0.002279036295,   0.000000000000,
-  -0.000743287324,   0.000000000000,
-   0.000567271195,   0.000000000000,
+  -0.000083604971,   0.000000000000,
+   0.001663378142,   0.000000000000,
+  -0.001020100652,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -528,31 +528,31 @@ const PylithScalar pylith::faults::CohesiveDynDataTri3d::_fieldIncrOpen[] = {
   1.2, 2.2, // 10
   1.3, 2.3, // 11
   1.5, 2.5, // 12
--10.8, 0.8, // 26
--10.0, 0.1, // 27
-  1.2,-10.2, // 28
++11.8, 11.8, // 26
++10.0, 0.1, // 27
+  1.2,+10.2, // 28
 };
 
 // Output
 const PylithScalar pylith::faults::CohesiveDynDataTri3d::_fieldIncrOpenE[] = {
    1.100000000000,   2.100000000000,
-   1.208740107962,   2.201462120349,
-   1.304388525058,   2.303014075601,
+   1.190062443638,   2.192265294477,
+   1.292981353233,   2.293412522606,
    1.400000000000,   2.400000000000,
-   1.502122067204,   2.503667076610,
+   1.495037703623,   2.494851226153,
    1.600000000000,   2.600000000000,
-   1.191259892038,   2.198537879651,
-   1.295611474942,   2.296985924399,
-   1.497877932796,   2.496332923390,
-   3.800000000000,  -4.800000000000,
-  -3.000000000000,  -4.000000000000,
-  -3.200000000000,  -4.200000000000,
+   1.209937556362,   2.207734705523,
+   1.307018646767,   2.306587477394,
+   1.504962296377,   2.505148773847,
+   3.800000000000,   4.800000000000,
+   3.000000000000,  -4.000000000000,
+  -3.200000000000,   4.200000000000,
 };
 
 const PylithScalar pylith::faults::CohesiveDynDataTri3d::_slipOpenE[] = {
-  -0.010292628789,   0.014428129643,
-   0.006028151203,   0.008777050116,
-  -0.004244134409,   0.007334153220,
+  -0.003115301533,   0.024992352435,
+   0.013174954788,   0.014037293534,
+  -0.009924592753,   0.010297547694,
 };
 
 // ----------------------------------------------------------------------
diff --git a/unittests/libtests/faults/data/cohesivedyn.py b/unittests/libtests/faults/data/cohesivedyn.py
index d493f7f..632ce98 100644
--- a/unittests/libtests/faults/data/cohesivedyn.py
+++ b/unittests/libtests/faults/data/cohesivedyn.py
@@ -1,5 +1,4 @@
-cell = "hex8"
-dim = "3d"
+cell = "tri3d"
 testCase = "open"
 
 import numpy
@@ -46,7 +45,7 @@ def faultToGlobal(v, R):
 
 
 # ----------------------------------------------------------------------
-if dim == "2d":
+if cell == "tri3" or cell == "tri3d" or cell == "quad4":
     if cell == "tri3":
         dlagrange1 = numpy.zeros(2)
         indexL = numpy.arange(12,16)
@@ -56,30 +55,30 @@ if dim == "2d":
         m = 4
         DOF = 2
 
-        fieldT = numpy.array([[8.6, 9.6],
-                              [8.8, 9.8]])
-        fieldIncr = numpy.array([[1.6, 2.6],
-                                 [1.8, 2.8]])
+        fieldT = numpy.array([[-8.6, 9.6],
+                              [-8.8, 9.8]])
+        fieldIncr = numpy.array([[-1.6, 2.6],
+                                 [-1.8, 2.8]])
         L = numpy.array([[1.0, 0.0, 0.0, 0.0,],
                          [0.0, 1.0, 0.0, 0.0,],
                          [0.0, 0.0, 1.0, 0.0,],
                          [0.0, 0.0, 0.0, 1.0,],]);
-        C = numpy.array([[0.0, -1.0, 0.0, 0.0,],
-                         [-1.0, 0.0, 0.0, 0.0,],
-                         [0.0, 0.0, 0.0, -1.0,],
-                         [0.0, 0.0, -1.0, 0.0,],]);
+        C = numpy.array([[0.0, +1.0, 0.0, 0.0,],
+                         [+1.0, 0.0, 0.0, 0.0,],
+                         [0.0, 0.0, 0.0, +1.0,],
+                         [0.0, 0.0, +1.0, 0.0,],]);
     
         jacobianN = numpy.array(
-            [[  4.0,  -1.2,  -2.2,  -2.3,],
-             [  -1.2,  5.0,  -1.3,  -3.2,],
-             [  -2.2,  -1.3,  4.1,  -4.3,],
-             [  -2.3,  -3.2,  -4.3,  5.1,],])
+            [[   4.0,  -1.2,  -2.2,  -2.3,],
+             [  -1.2,   5.0,  -1.3,  -3.2,],
+             [  -2.2,  -1.3,   4.1,  -4.3,],
+             [  -2.3,  -3.2,  -4.3,   5.1,],])
 
         jacobianP = numpy.array(
-            [[  5.0,  -1.2,  -2.2,  -2.3,],
-             [  -1.2,  4.0,  -1.3,  -3.2,],
-             [  -2.2,  -1.3,  5.1,  -4.3,],
-             [  -2.3,  -3.2,  -4.3,  4.1,],])
+            [[   5.0,  -1.2,  -2.2,  -2.3,],
+             [  -1.2,   4.0,  -1.3,  -3.2,],
+             [  -2.2,  -1.3,   5.1,  -4.3,],
+             [  -2.3,  -3.2,  -4.3,   4.1,],])
 
         disp = numpy.array([[ 8.1, 9.1,],
                             [ 8.2, 9.2,],
@@ -87,8 +86,8 @@ if dim == "2d":
                             [ 8.4, 9.4,],
                             [ 8.2, 9.2,],
                             [ 8.3, 9.3,],
-                            [ 8.6, 9.6,],
-                            [ 8.8, 9.8,],])
+                            [-8.6, 9.6,],
+                            [-8.8, 9.8,],])
 
         if testCase == "slip":
             dispIncr = numpy.array([[ 9.1, 7.1,],
@@ -97,8 +96,8 @@ if dim == "2d":
                                     [ 9.4, 7.4,],
                                     [ 9.2, 7.2,],
                                     [ 9.3, 7.3,],
-                                    [ 1.6, 2.6,],
-                                    [ 1.8, 2.8,],])            
+                                    [-1.6, 2.6,],
+                                    [-1.8, 2.8,],])            
         elif testCase == "open":
             dispIncr = numpy.array([[ 9.1, 7.1,],
                                     [ 9.2, 7.2,],
@@ -106,8 +105,8 @@ if dim == "2d":
                                     [ 9.4, 7.4,],
                                     [ 9.2, 7.2,],
                                     [ 9.3, 7.3,],
-                                    [ -10.6, 2.6,],
-                                    [ -10.8, 2.8,],])            
+                                    [ +10.6, -10.6,],
+                                    [ +10.8, -10.8,],])            
 
 
     elif cell == "tri3d":
@@ -119,12 +118,12 @@ if dim == "2d":
         m = 6
         DOF = 2
 
-        fieldT = numpy.array([[-3.8, 4.8],
-                              [3.0, 4.0],
-                              [3.2, 4.2]])
-        fieldIncr = numpy.array([[1.8, 0.8],
-                                 [1.0, 0.1],
-                                 [1.2, 0.2]])
+        fieldT = numpy.array([[-3.8,-4.8],
+                              [-3.0, 4.0],
+                              [3.2, -4.2]])
+        fieldIncr = numpy.array([[-1.8,+3.6],
+                                 [-1.0, 1.1],
+                                 [ 1.7,-1.2]])
 
         L = numpy.array([[2.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                          [0.0, 2.0, 0.0, 0.0, 0.0, 0.0],
@@ -133,12 +132,12 @@ if dim == "2d":
                          [0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
                          [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]])
 
-        C = numpy.array([[+0.70710678118654757, -0.70710678118654757, 0.0, 0.0, 0.0, 0.0,],
-                         [-0.70710678118654757, -0.70710678118654757, 0.0, 0.0, 0.0, 0.0,],
-                         [0.0, 0.0, 0.0, -1.0, 0.0, 0.0,],
-                         [0.0, 0.0, -1.0, 0.0, 0.0, 0.0,],
-                         [0.0, 0.0, 0.0, 0.0, +1.0, 0.0,],
-                         [0.0, 0.0, 0.0, 0.0, 0.0, -1.0,],])
+        C = numpy.array([[-0.70710678118654757, +0.70710678118654757, 0.0, 0.0, 0.0, 0.0,],
+                         [+0.70710678118654757, +0.70710678118654757, 0.0, 0.0, 0.0, 0.0,],
+                         [0.0, 0.0, 0.0, +1.0, 0.0, 0.0,],
+                         [0.0, 0.0, +1.0, 0.0, 0.0, 0.0,],
+                         [0.0, 0.0, 0.0, 0.0, -1.0, 0.0,],
+                         [0.0, 0.0, 0.0, 0.0, 0.0, +1.0,],])
     
         jacobianN = numpy.array(
             [[+6.0, -1.0, -1.1, -1.2, -1.3, -1.4],
@@ -165,9 +164,9 @@ if dim == "2d":
                             [ 6.2, 8.2,],
                             [ 6.3, 8.3,],
                             [ 6.5, 8.5,],
-                            [-3.8, 4.8,],
-                            [ 3.0, 4.0,],
-                            [ 3.2, 4.2,],])
+                            [-3.8,-4.8,],
+                            [-3.0, 4.0,],
+                            [ 3.2,-4.2,],])
 
         if testCase == "slip":
             dispIncr = numpy.array([[ 1.1, 2.1,],
@@ -179,9 +178,9 @@ if dim == "2d":
                                     [ 1.2, 2.2,],
                                     [ 1.3, 2.3,],
                                     [ 1.5, 2.5,],
-                                    [ 1.8, 0.8,],
-                                    [ 1.0, 0.1,],
-                                    [ 1.2, 0.2,],])            
+                                    [-1.8,+3.6,],
+                                    [-1.0, 1.1,],
+                                    [ 1.7,-1.2,],])            
         elif testCase == "open":
             dispIncr = numpy.array([[ 1.1, 2.1,],
                                     [ 1.2, 2.2,],
@@ -192,9 +191,9 @@ if dim == "2d":
                                     [ 1.2, 2.2,],
                                     [ 1.3, 2.3,],
                                     [ 1.5, 2.5,],
-                                    [-10.8, 0.8,],
-                                    [-10.0, 0.1,],
-                                    [ 1.2, -10.2,],])            
+                                    [+11.8, 11.8,],
+                                    [+10.0, 0.1,],
+                                    [ 1.2, +10.2,],])            
 
 
     elif cell == "quad4":
@@ -206,64 +205,64 @@ if dim == "2d":
         m = 4
         DOF = 2
 
-        fieldT = numpy.array([[8.8, 9.8],
-                              [8.0, 9.0]])
-        fieldIncr = numpy.array([[1.8, 2.8],
-                                 [1.0, 2.0]])
+        fieldT = numpy.array([[-8.6, 9.6],
+                              [-8.8, 9.8]])
+        fieldIncr = numpy.array([[-1.6, 2.6],
+                                 [-1.8, 2.5]])
         L = numpy.array([[1.0, 0.0, 0.0, 0.0,],
                          [0.0, 1.0, 0.0, 0.0,],
                          [0.0, 0.0, 1.0, 0.0,],
                          [0.0, 0.0, 0.0, 1.0,],]);
-        C = numpy.array([[0.0, -1.0, 0.0, 0.0,],
-                         [-1.0, 0.0, 0.0, 0.0,],
-                         [0.0, 0.0, 0.0, -1.0,],
-                         [0.0, 0.0, -1.0, 0.0,],]);
+        C = numpy.array([[0.0, +1.0, 0.0, 0.0,],
+                         [+1.0, 0.0, 0.0, 0.0,],
+                         [0.0, 0.0, 0.0, +1.0,],
+                         [0.0, 0.0, +1.0, 0.0,],]);
     
         jacobianN = numpy.array(
-            [[  4.0,  -1.2,  -2.2,  -2.3,],
-             [  -1.2,  5.0,  -1.3,  -3.2,],
-             [  -2.2,  -1.3,  4.1,  -4.3,],
-             [  -2.3,  -3.2,  -4.3,  5.1,],])
+            [[   4.0,  -1.2,  -2.2,  -2.3,],
+             [  -1.2,   5.0,  -1.3,  -3.2,],
+             [  -2.2,  -1.3,   4.1,  -4.3,],
+             [  -2.3,  -3.2,  -4.3,   5.1,],])
 
         jacobianP = numpy.array(
-            [[  5.0,  -1.2,  -2.2,  -2.3,],
-             [  -1.2,  4.0,  -1.3,  -3.2,],
-             [  -2.2,  -1.3,  5.1,  -4.3,],
-             [  -2.3,  -3.2,  -4.3,  4.1,],])
+            [[   5.0,  -1.2,  -2.2,  -2.3,],
+             [  -1.2,   4.0,  -1.3,  -3.2,],
+             [  -2.2,  -1.3,   5.1,  -4.3,],
+             [  -2.3,  -3.2,  -4.3,   4.1,],])
 
         disp = numpy.array([[ 8.1, 9.1,],
+                            [ 8.3, 9.3,],
                             [ 8.2, 9.2,],
                             [ 8.3, 9.3,],
-                            [ 8.4, 9.4,],
                             [ 8.5, 9.5,],
                             [ 8.6, 9.6,],
-                            [ 8.3, 9.3,],
-                            [ 8.4, 9.4,],
-                            [ 8.8, 9.8,],
-                            [ 8.0, 9.0,],])
+                            [ 8.2, 9.6,],
+                            [ 8.3, 9.8,],
+                            [-8.6, 9.6,],
+                            [-8.8, 9.8,],])
 
         if testCase == "slip":
             dispIncr = numpy.array([[ 1.1, 2.1,],
                                     [ 1.2, 2.2,],
+                                    [ 1.2, 2.2,],
                                     [ 1.3, 2.3,],
-                                    [ 1.4, 2.4,],
                                     [ 1.5, 2.5,],
                                     [ 1.6, 2.6,],
+                                    [ 1.2, 2.2,],
                                     [ 1.3, 2.3,],
-                                    [ 1.4, 2.4,],
-                                    [ 1.8, 2.8,],
-                                    [ 1.0, 2.0,],])            
+                                    [-1.6, 2.6,],
+                                    [-1.8, 2.5,],])            
         elif testCase == "open":
             dispIncr = numpy.array([[ 1.1, 2.1,],
                                     [ 1.2, 2.2,],
+                                    [ 1.2, 2.2,],
                                     [ 1.3, 2.3,],
-                                    [ 1.4, 2.4,],
                                     [ 1.5, 2.5,],
                                     [ 1.6, 2.6,],
+                                    [ 1.2, 2.2,],
                                     [ 1.3, 2.3,],
-                                    [ 1.4, 2.4,],
-                                    [ -10.8, 2.8,],
-                                    [ -10.0, 2.0,],])            
+                                    [+10.6, -12.6,],
+                                    [+10.8, -12.8,],])            
 
 
     # ------------------------------------------------------------------
@@ -303,21 +302,21 @@ if dim == "2d":
     duN = numpy.dot(inv(jacobianN),RHS)
     duP = -numpy.dot(inv(jacobianP),RHS)
     
-    dispRel = duP - duN
+    dispRelIncr = duP - duN
 
     dispTpdt = disp + dispIncr
     dispTpdt = numpy.reshape(dispTpdt, n)
 
-    slipVertex = dispRel + dispTpdt[indexP]-dispTpdt[indexN]
+    slipVertex = dispRelIncr + dispTpdt[indexP]-dispTpdt[indexN]
     slipVertex = numpy.reshape(slipVertex, (m/DOF,DOF))
     slipVertex = globalToFault(slipVertex, C)
-    if testCase == "slip":
-        slipVertex[:,1] = 0
     mask = slipVertex[:,1] < 0.0
-    slipVertex[mask,1] = 0
+    #slipVertex[:,1] = 0
     print "slip",slipVertex
     slipVertex = faultToGlobal(slipVertex, C)
     slipVertex = numpy.reshape(slipVertex, m)
+    disp = numpy.reshape(disp, n)
+    slipIncrVertex = slipVertex - (disp[indexP] - disp[indexN])
 
     print "duN \n", duN
     print "duP \n", duP
@@ -325,8 +324,8 @@ if dim == "2d":
     dispIncrE = dispIncr
     dispIncrE = numpy.reshape(dispIncrE, n)
     dispIncrE[indexL] = dispIncrE[indexL] + dLagrange
-    dispIncrE[indexN] = dispIncrE[indexN] - 0.5*slipVertex
-    dispIncrE[indexP] = dispIncrE[indexP] + 0.5*slipVertex
+    dispIncrE[indexN] = dispIncrE[indexN] - 0.5*slipIncrVertex
+    dispIncrE[indexP] = dispIncrE[indexP] + 0.5*slipIncrVertex
     dispIncrE = numpy.reshape(dispIncrE, (n/DOF,DOF))
 
     slipVertex = numpy.reshape(slipVertex, (m/DOF,DOF))
@@ -337,7 +336,7 @@ if dim == "2d":
 
 
 # ----------------------------------------------------------------------
-elif dim == "3d":
+elif cell == "tet4" or cell == "hex8":
     lagrangeScale = lengthScale**2
 
     if cell == "tet4":
@@ -350,12 +349,12 @@ elif dim == "3d":
         m = 9
         DOF = 3
 
-        fieldT = numpy.array([[7.7, 8.7, 9.7],
-                              [7.9, 8.9, 9.9],
-                              [7.1, 8.1, 9.1]])
-        fieldIncr = numpy.array([[9.7, 2.7, 3.7],
-                                 [9.9, 2.9, 3.9],
-                                 [9.1, 2.1, 3.1]])
+        fieldT = numpy.array([[-7.7, 18.7, 19.7],
+                              [-7.9, 18.9, 19.9],
+                              [-7.1, 18.1, 19.1]])
+        fieldIncr = numpy.array([[-4.7, 5.7, 6.7],
+                                 [-4.9, 5.9, 6.9],
+                                 [-4.1, 5.1, 6.1]])
         
         L = numpy.array([[1.0/3.0,0,0, 0.0,0,0, 0.0,0,0,],
                          [0,1.0/3.0,0, 0,0.0,0, 0,0.0,0,],
@@ -367,9 +366,9 @@ elif dim == "3d":
                          [0,0.0,0, 0,0.0,0, 0,1.0/3.0,0,],
                          [0,0,0.0, 0,0,0.0, 0,0,1.0/3.0,]])
 
-        Cv = numpy.array([[ 0, -1, 0,],
+        Cv = numpy.array([[ 0, +1, 0,],
                           [ 0, 0, +1,],
-                          [ -1, 0, 0,],])
+                          [ +1, 0, 0,],])
         Zv = numpy.zeros([3,3])
         C = numpy.vstack( (numpy.hstack((Cv, Zv, Zv)),
                            numpy.hstack((Zv, Cv, Zv)),
@@ -405,9 +404,9 @@ elif dim == "3d":
                             [ 7.2, 8.2, 9.2,],
                             [ 7.3, 8.3, 9.3,],
                             [ 7.4, 8.4, 9.4,],
-                            [ 7.7, 8.7, 9.7,],
-                            [ 7.9, 8.9, 9.9,],
-                            [ 7.1, 8.1, 9.1,],])
+                            [-7.7, 18.7, 19.7,],
+                            [-7.9, 18.9, 19.9,],
+                            [-7.1, 18.1, 19.1,],])
 
         if testCase == "slip":
             dispIncr = numpy.array([[ 1.1, 2.1, 3.1,],
@@ -418,9 +417,9 @@ elif dim == "3d":
                                     [ 1.2, 2.2, 3.2,],
                                     [ 1.3, 2.3, 3.3,],
                                     [ 1.4, 2.4, 3.4,],
-                                    [ 9.7, 2.7, 3.7,],
-                                    [ 9.9, 2.9, 3.9,],
-                                    [ 9.1, 2.1, 3.1,],])            
+                                    [-4.7, 5.7, 6.7,],
+                                    [-4.9, 5.9, 6.9,],
+                                    [-4.1, 5.1, 6.1,],])            
         elif testCase == "open":
             dispIncr = numpy.array([[ 1.1, 2.1, 3.1,],
                                     [ 1.2, 2.2, 3.2,],
@@ -430,9 +429,9 @@ elif dim == "3d":
                                     [ 1.2, 2.2, 3.2,],
                                     [ 1.3, 2.3, 3.3,],
                                     [ 1.4, 2.4, 3.4,],
-                                    [-20.7,  2.7, 3.7,],
-                                    [-20.9,  2.9, 3.9,],
-                                    [-20.1,  2.1, 3.1,],])            
+                                    [+80.7,  2.7, 3.7,],
+                                    [+80.9,  2.9, 3.9,],
+                                    [+80.1,  2.1, 3.1,],])            
 
 
     elif cell == "hex8":
@@ -460,18 +459,18 @@ elif dim == "3d":
                          [0, a2, 0, 0, a1, 0, 0, a1, 0, 0, a0, 0],
                          [0, 0, a2, 0, 0, a1, 0, 0, a1, 0, 0, a0]])
 
-        fieldT = numpy.array([[4.4, 2.4, 3.4],
-                              [4.6, 2.6, 3.6],
-                              [4.8, 2.8, 3.8],
-                              [4.0, 2.0, 3.0]])
-        fieldIncr = numpy.array([[1.4, 2.4, 0.4],
-                                 [1.6, 2.6, 0.6],
-                                 [1.8, 2.8, 0.8],
-                                 [1.0, 2.0, 0.2]])
+        fieldT = numpy.array([[-4.4, 2.4, 3.4],
+                              [-4.6, 2.6, 3.6],
+                              [-4.8, 2.8, 3.8],
+                              [-4.0, 2.0, 3.0]])
+        fieldIncr = numpy.array([[-1.4, 2.4, 0.4],
+                                 [-1.6, 2.6, 0.6],
+                                 [-1.8, 2.8, 0.8],
+                                 [-1.0, 2.0, 0.2]])
 
-        Cv = numpy.array([[ 0, -1, 0,],
+        Cv = numpy.array([[ 0, +1, 0,],
                           [ 0, 0, +1,],
-                          [ -1, 0, 0,],])
+                          [ +1, 0, 0,],])
         Zv = numpy.zeros([3,3])
         C = numpy.vstack( (numpy.hstack((Cv, Zv, Zv, Zv)),
                            numpy.hstack((Zv, Cv, Zv, Zv)),
@@ -519,14 +518,14 @@ elif dim == "3d":
                             [ 4.0, 2.0, 3.0,],
                             [ 4.1, 2.1, 3.1,],
                             [ 4.2, 2.2, 3.2,],
-                            [ 4.5, 2.5, 3.5,],
-                            [ 4.6, 2.6, 3.6,],
-                            [ 4.7, 2.7, 3.7,],
-                            [ 4.8, 2.8, 3.8,],
-                            [ 4.4, 2.4, 3.4,],
-                            [ 4.6, 2.6, 3.6,],
-                            [ 4.8, 2.8, 3.8,],
-                            [ 4.0, 2.0, 3.0,],])
+                            [ 4.5, 3.2, 4.3,],
+                            [ 4.6, 3.5, 4.6,],
+                            [ 4.7, 3.7, 4.6,],
+                            [ 4.8, 3.6, 4.5,],
+                            [-4.4, 2.4, 3.4,],
+                            [-4.6, 2.6, 3.6,],
+                            [-4.8, 2.8, 3.8,],
+                            [-4.0, 2.0, 3.0,],])
 
         if testCase == "slip":
             dispIncr = numpy.array([[ 1.1, 2.1, 0.1,],
@@ -545,10 +544,11 @@ elif dim == "3d":
                                     [ 1.6, 2.6, 0.6,],
                                     [ 1.7, 2.7, 0.7,],
                                     [ 1.8, 2.8, 0.8,],
-                                    [ 1.4, 2.4, 0.4,],
-                                    [ 1.6, 2.6, 0.6,],
-                                    [ 1.8, 2.8, 0.8,],
-                                    [ 1.0, 2.0, 0.2,],])          
+                                    [-1.4, 2.4, 0.4,],
+                                    [-1.6, 2.6, 0.6,],
+                                    [-1.8, 2.8, 0.8,],
+                                    [-1.0, 2.0, 0.2,],
+                                    ])
         elif testCase == "open":
             dispIncr = numpy.array([[ 1.1, 2.1, 0.1,],
                                     [ 1.2, 2.2, 0.2,],
@@ -566,10 +566,10 @@ elif dim == "3d":
                                     [ 1.6, 2.6, 0.6,],
                                     [ 1.7, 2.7, 0.7,],
                                     [ 1.8, 2.8, 0.8,],
-                                    [-10.4, 2.4, 0.4,],
-                                    [-10.6, 2.6, 0.6,],
-                                    [-10.8, 2.8, 0.8,],
-                                    [-10.0, 2.0, 0.2,],])          
+                                    [+20.4, 2.4, 0.4,],
+                                    [+20.6, 2.6, 0.6,],
+                                    [+20.8, 2.8, 0.8,],
+                                    [+20.0, 2.0, 0.2,],])          
 
     # ------------------------------------------------------------------
     fieldTpdt = fieldT + fieldIncr
@@ -622,6 +622,8 @@ elif dim == "3d":
     slipVertex[mask,2] = 0
     slipVertex = faultToGlobal(slipVertex, C)
     slipVertex = numpy.reshape(slipVertex, m)
+    disp = numpy.reshape(disp, n)
+    slipIncrVertex = slipVertex - (disp[indexP] - disp[indexN])
 
     print "duN \n", duN
     print "duP \n", duP
@@ -630,8 +632,8 @@ elif dim == "3d":
     dispIncrE = numpy.reshape(dispIncrE, n)
 
     dispIncrE[indexL] = dispIncrE[indexL] + dLagrange
-    dispIncrE[indexN] = dispIncrE[indexN] - 0.5*slipVertex
-    dispIncrE[indexP] = dispIncrE[indexP] + 0.5*slipVertex
+    dispIncrE[indexN] = dispIncrE[indexN] - 0.5*slipIncrVertex
+    dispIncrE[indexP] = dispIncrE[indexP] + 0.5*slipIncrVertex
     dispIncrE = numpy.reshape(dispIncrE, (n/DOF,DOF))
 
     slipVertex = numpy.reshape(slipVertex, (m/DOF,DOF))



More information about the CIG-COMMITS mailing list