[cig-commits] [commit] baagaard/dynrup-new-lagrange: More debugging of new FaultCohesiveDyn using unit tests. (88ac92e)

cig_noreply at geodynamics.org cig_noreply at geodynamics.org
Fri Jun 6 17:32:23 PDT 2014


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

On branch  : baagaard/dynrup-new-lagrange
Link       : https://github.com/geodynamics/pylith/compare/ba6825e04f1851181ee6fadc84526954d3e896b8...88ac92e5aca24f5a10bce154bb9b865699eed75a

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

commit 88ac92e5aca24f5a10bce154bb9b865699eed75a
Author: Brad Aagaard <baagaard at usgs.gov>
Date:   Fri Jun 6 17:32:14 2014 -0700

    More debugging of new FaultCohesiveDyn using unit tests.
    
    Started updating tri3 test data using cohesivedyn.py.


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

88ac92e5aca24f5a10bce154bb9b865699eed75a
 unittests/libtests/faults/TestFaultCohesiveDyn.cc  |  79 ++++----
 .../libtests/faults/data/CohesiveDynDataTri3.cc    |  75 ++++----
 unittests/libtests/faults/data/cohesivedyn.py      | 199 +++++++--------------
 3 files changed, 141 insertions(+), 212 deletions(-)

diff --git a/unittests/libtests/faults/TestFaultCohesiveDyn.cc b/unittests/libtests/faults/TestFaultCohesiveDyn.cc
index 81418c3..a5e2880 100644
--- a/unittests/libtests/faults/TestFaultCohesiveDyn.cc
+++ b/unittests/libtests/faults/TestFaultCohesiveDyn.cc
@@ -243,15 +243,15 @@ pylith::faults::TestFaultCohesiveDyn::testIntegrateResidualStick(void)
   residual.view("RESIDUAL"); // DEBUGGING
   { // Check residual values
     PetscInt pStart=0, pEnd=0;
-    topology::VecVisitorMesh residualVisitor(residual);
     PetscErrorCode err = PetscSectionGetChart(residual.petscSection(), &pStart, &pEnd);CPPUNIT_ASSERT(!err);
+    topology::VecVisitorMesh residualVisitor(residual);
     const PetscScalar* residualArray = residualVisitor.localArray();CPPUNIT_ASSERT(residualArray);
 
     const PylithScalar* residualE = _data->residualStickE;CPPUNIT_ASSERT(residualE);
     const int fiberDimE = spaceDim; // number of values per point
     const PylithScalar tolerance = 1.0e-06;
 
-    for (PetscInt p=pStart, iPoint=0; p < pEnd; ++p, ++iPoint) {
+    for (PetscInt p=pStart, iPoint=0; p < pEnd; ++p) {
       if (residualVisitor.sectionDof(p) > 0) {
 	const PetscInt off = residualVisitor.sectionOffset(p);
 	CPPUNIT_ASSERT_EQUAL(fiberDimE, residualVisitor.sectionDof(p));
@@ -263,6 +263,7 @@ pylith::faults::TestFaultCohesiveDyn::testIntegrateResidualStick(void)
 	    CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, residualArray[off+iDim], tolerance);
 	  } // if/else
 	} // for
+	++iPoint;
       } // if
     } // for
   } // Check residual values
@@ -295,30 +296,31 @@ pylith::faults::TestFaultCohesiveDyn::testIntegrateResidualSlip(void)
   const topology::Field& residual = fields.get("residual");
   fault.integrateResidual(residual, t, &fields);
 
-  //residual.view("RESIDUAL"); // DEBUGGING
+  residual.view("RESIDUAL"); // DEBUGGING
   { // Check residual values
-    PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
-    topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
-    const PetscInt vStart = verticesStratum.begin();
-    const PetscInt vEnd = verticesStratum.end();
-
-    topology::VecVisitorMesh residualVisitor(fields.get("residual"));
+    PetscInt pStart=0, pEnd=0;
+    PetscErrorCode err = PetscSectionGetChart(residual.petscSection(), &pStart, &pEnd);CPPUNIT_ASSERT(!err);
+    topology::VecVisitorMesh residualVisitor(residual);
     const PetscScalar* residualArray = residualVisitor.localArray();CPPUNIT_ASSERT(residualArray);
 
     const PylithScalar* residualE = _data->residualSlipE;CPPUNIT_ASSERT(residualE);
     const int fiberDimE = spaceDim; // number of values per point
     const PylithScalar tolerance = 1.0e-06;
-    for(PetscInt v = vStart, iVertex = 0; v < vEnd; ++v, ++iVertex) {
-      const PetscInt off = residualVisitor.sectionOffset(v);
-      CPPUNIT_ASSERT_EQUAL(fiberDimE, residualVisitor.sectionDof(v));
-      for(PetscInt d = 0; d < fiberDimE; ++d) {
-        const PylithScalar valE = residualE[iVertex*spaceDim+d];
-        if (fabs(valE) > tolerance) {
-          CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, residualArray[off+d]/valE, tolerance);
-        } else {
-          CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, residualArray[off+d], tolerance);
-	} // if/else
-      } // for
+
+    for (PetscInt p=pStart, iPoint=0; p < pEnd; ++p) {
+      if (residualVisitor.sectionDof(p) > 0) {
+	const PetscInt off = residualVisitor.sectionOffset(p);
+	CPPUNIT_ASSERT_EQUAL(fiberDimE, residualVisitor.sectionDof(p));
+	for (PetscInt iDim=0; iDim < fiberDimE; ++iDim) {
+	  const PylithScalar valE = residualE[iPoint*fiberDimE+iDim];
+	  if (fabs(valE) > tolerance) {
+	    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, residualArray[off+iDim]/valE, tolerance);
+	  } else {
+	    CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, residualArray[off+iDim], tolerance);
+	  } // if/else
+	} // for
+	++iPoint;
+      } // if
     } // for
   } // Check residual values
 
@@ -356,30 +358,31 @@ pylith::faults::TestFaultCohesiveDyn::testIntegrateResidualOpen(void)
 
   fault.updateStateVars(t, &fields);
 
-  //residual.view("RESIDUAL"); // DEBUGGING
+  residual.view("RESIDUAL"); // DEBUGGING
   { // Check residual values
-    PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
-    topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
-    const PetscInt vStart = verticesStratum.begin();
-    const PetscInt vEnd = verticesStratum.end();
-
-    topology::VecVisitorMesh residualVisitor(fields.get("residual"));
+    PetscInt pStart=0, pEnd=0;
+    PetscErrorCode err = PetscSectionGetChart(residual.petscSection(), &pStart, &pEnd);CPPUNIT_ASSERT(!err);
+    topology::VecVisitorMesh residualVisitor(residual);
     const PetscScalar* residualArray = residualVisitor.localArray();CPPUNIT_ASSERT(residualArray);
 
     const PylithScalar* residualE = _data->residualOpenE;CPPUNIT_ASSERT(residualE);
     const int fiberDimE = spaceDim; // number of values per point
     const PylithScalar tolerance = 1.0e-06;
-    for(PetscInt v = vStart, iVertex = 0; v < vEnd; ++v, ++iVertex) {
-      const PetscInt off = residualVisitor.sectionOffset(v);
-      CPPUNIT_ASSERT_EQUAL(fiberDimE, residualVisitor.sectionDof(v));
-      for(PetscInt d = 0; d < fiberDimE; ++d) {
-        const PylithScalar valE = residualE[iVertex*spaceDim+d];
-        if (fabs(valE) > tolerance) {
-          CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, residualArray[off+d]/valE, tolerance);
-        } else {
-          CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, residualArray[off+d], tolerance);
-	} // if/else
-      } // for
+
+    for (PetscInt p=pStart, iPoint=0; p < pEnd; ++p) {
+      if (residualVisitor.sectionDof(p) > 0) {
+	const PetscInt off = residualVisitor.sectionOffset(p);
+	CPPUNIT_ASSERT_EQUAL(fiberDimE, residualVisitor.sectionDof(p));
+	for (PetscInt iDim=0; iDim < fiberDimE; ++iDim) {
+	  const PylithScalar valE = residualE[iPoint*fiberDimE+iDim];
+	  if (fabs(valE) > tolerance) {
+	    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, residualArray[off+iDim]/valE, tolerance);
+	  } else {
+	    CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, residualArray[off+iDim], tolerance);
+	  } // if/else
+	} // for
+	++iPoint;
+      } // if
     } // for
   } // Check residual values
 
diff --git a/unittests/libtests/faults/data/CohesiveDynDataTri3.cc b/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
index 5f1e7de..8ec71ce 100644
--- a/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
+++ b/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
@@ -240,17 +240,17 @@ const PylithScalar pylith::faults::CohesiveDynDataTri3::_jacobian[] = {
 // Computed values
 // ----------------------------------------------------------------------
 
-const PylithScalar pylith::faults::CohesiveDynDataTri3::_orientation[] = {
+const PylithScalar pylith::faults::CohesiveDynDataTri3::_orientation[2*2*2] = {
   0.0, +1.0,  +1.0, 0.0,
   0.0, +1.0,  +1.0, 0.0
 };
 
-const PylithScalar pylith::faults::CohesiveDynDataTri3::_area[] = {
+const PylithScalar pylith::faults::CohesiveDynDataTri3::_area[2] = {
   1.0,
   1.0,
 };
 
-const PylithScalar pylith::faults::CohesiveDynDataTri3::_initialTractions[] = {
+const PylithScalar pylith::faults::CohesiveDynDataTri3::_initialTractions[2*2] = {
   // Fault coordinate frame
   1.0, -2.0,
   1.1, -2.1,
@@ -258,18 +258,19 @@ const PylithScalar pylith::faults::CohesiveDynDataTri3::_initialTractions[] = {
 
 
 const int pylith::faults::CohesiveDynDataTri3::_numConstraintEdges = 2;
-const int pylith::faults::CohesiveDynDataTri3::_constraintEdges[] = {
+const int pylith::faults::CohesiveDynDataTri3::_constraintEdges[2] = {
   15, 16
 };
-const int pylith::faults::CohesiveDynDataTri3::_negativeVertices[] = {
+const int pylith::faults::CohesiveDynDataTri3::_negativeVertices[2] = {
   4, 5
 };
 
 // ----------------------------------------------------------------------
 // Stick case
 // ----------------------------------------------------------------------
+// Computed using cohesivedyn.py
 // Input
-const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrStick[] = {
+const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrStick[8*2] = {
   1.1, 2.1,
   1.2, 2.2, // 3
   1.3, 2.3, // 4
@@ -281,42 +282,42 @@ const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrStick[] = {
 };
 
 // Output
-const PylithScalar pylith::faults::CohesiveDynDataTri3::_residualStickE[] = {
-  0.0, 0.0,
-  0.0, 0.001*(0.6*(21.6+8.6)-(9.6+2.6)), // 3
-  0.0, 0.001*(0.6*(21.8+8.8)-(9.8+2.8)), // 4
-  0.0, 0.0,
-  0.0, -0.001*(0.6*(21.6+8.6)-(9.6+2.6)), // 6
-  0.0, -0.001*(0.6*(21.8+8.8)-(9.8+2.8)), // 7
- -21.6, 2.6, // 15
- -21.8, 2.8, // 16
+const PylithScalar pylith::faults::CohesiveDynDataTri3::_residualStickE[8*2] = {
+   0.0000000000e+00,   0.0000000000e+00,
+   0.0000000000e+00,   5.9200000000e-03,
+   0.0000000000e+00,   5.7600000000e-03,
+   0.0000000000e+00,   0.0000000000e+00,
+  -0.0000000000e+00,  -5.9200000000e-03,
+  -0.0000000000e+00,  -5.7600000000e-03,
+  -0.0000000000e+00,   0.0000000000e+00,
+  -0.0000000000e+00,   0.0000000000e+00,
 };
 
 // ----------------------------------------------------------------------
 // Slip case
 // ----------------------------------------------------------------------
 // Input
-const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrSlip[] = {
+const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrSlip[8*2] = {
   9.1, 7.1,
   9.2, 7.2, // 3
   9.3, 7.3, // 4
   9.4, 7.4,
-  9.2, 7.2, // 6
-  9.3, 7.3, // 7
+  9.2, 7.5, // 6
+  9.3, 7.6, // 7
  -1.6, 2.6, // 15
  -1.8, 2.8, // 16
 };
 
 // Output
-const PylithScalar pylith::faults::CohesiveDynDataTri3::_residualSlipE[] = {
-   9.100000000000,   7.100000000000,
-   9.200000000000,   7.200190546440,
-   9.300000000000,   7.300983993112,
-   9.400000000000,   7.400000000000,
-   9.200000000000,   7.199809453560,
-   9.300000000000,   7.299016006888,
-  -1.600000000000,  -3.480000000000,
-  -1.800000000000,  -3.440000000000,
+const PylithScalar pylith::faults::CohesiveDynDataTri3::_residualSlipE[8*2] = {
+   0.0000000000e+00,   0.0000000000e+00,
+   0.0000000000e+00,  -6.0800000000e-03,
+   0.0000000000e+00,  -6.2400000000e-03,
+   0.0000000000e+00,   0.0000000000e+00,
+  -0.0000000000e+00,   6.0800000000e-03,
+  -0.0000000000e+00,   6.2400000000e-03,
+  -0.0000000000e+00,   3.6600000000e-03,
+  -0.0000000000e+00,   3.7800000000e-03,
 };
 
 // ----------------------------------------------------------------------
@@ -328,22 +329,22 @@ const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrOpen[] = {
   9.2, 7.2, // 3
   9.3, 7.3, // 4
   9.4, 7.4,
-  9.2, 7.2, // 6
-  9.3, 7.3, // 7
+  9.5, 7.5, // 6
+  9.6, 7.6, // 7
   +10.6, -10.6, // 15
   +10.8, -10.8, // 16
 };
 
 // Output
 const PylithScalar pylith::faults::CohesiveDynDataTri3::_residualOpenE[] = {
-   9.100000000000,   7.100000000000,
-   9.199826714877,   7.200414257705,
-   9.299943020362,   7.300360547702,
-   9.400000000000,   7.400000000000,
-   9.200173285123,   7.199585742295,
-   9.300056979638,   7.299639452298,
-   8.600000000000,  -9.600000000000,
-   8.800000000000,  -9.800000000000,
+   0.0000000000e+00,   0.0000000000e+00,
+   0.0000000000e+00,  -2.0000000000e-04,
+   0.0000000000e+00,  -2.0000000000e-04,
+   0.0000000000e+00,   0.0000000000e+00,
+  -0.0000000000e+00,   2.0000000000e-04,
+  -0.0000000000e+00,   2.0000000000e-04,
+   6.0000000000e-04,  -3.0000000000e-04,
+   6.0000000000e-04,  -3.0000000000e-04,
 };
 
 // ----------------------------------------------------------------------
diff --git a/unittests/libtests/faults/data/cohesivedyn.py b/unittests/libtests/faults/data/cohesivedyn.py
index 632ce98..8fa3892 100644
--- a/unittests/libtests/faults/data/cohesivedyn.py
+++ b/unittests/libtests/faults/data/cohesivedyn.py
@@ -1,5 +1,5 @@
-cell = "tri3d"
-testCase = "open"
+cell = "tri3"
+testCase = "slip"
 
 import numpy
 
@@ -16,7 +16,7 @@ def printdata(data):
     Print data as C array.
     """
     (nrows, ncols) = data.shape
-    style = " %16.12f,"*ncols
+    style = " %18.10e,"*ncols
     for row in xrange(nrows):
         print (style % tuple(data[row,:]))
     return
@@ -29,7 +29,7 @@ def globalToFault(v, R):
     """
     (m,ndof) = v.shape
 
-    vF = numpy.dot(C, v.reshape(m*ndof,1))
+    vF = numpy.dot(R, v.reshape(m*ndof,1))
     return vF.reshape((m, ndof))
 
 
@@ -40,25 +40,20 @@ def faultToGlobal(v, R):
     """
     (m,ndof) = v.shape
 
-    vG = numpy.dot(C.transpose(), v.reshape(m*ndof,1))
+    vG = numpy.dot(R.transpose(), v.reshape(m*ndof,1))
     return vG.reshape((m, ndof))
 
 
 # ----------------------------------------------------------------------
 if cell == "tri3" or cell == "tri3d" or cell == "quad4":
     if cell == "tri3":
-        dlagrange1 = numpy.zeros(2)
-        indexL = numpy.arange(12,16)
-        indexN = numpy.arange(2,6)
-        indexP = numpy.arange(8,12)
-        n = 16
-        m = 4
+        indexL = numpy.arange(6,8)
+        indexN = numpy.arange(1,3)
+        indexP = numpy.arange(4,6)
+        n = 8
+        m = 2
         DOF = 2
 
-        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,],
@@ -68,45 +63,42 @@ if cell == "tri3" or cell == "tri3d" or cell == "quad4":
                          [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,],])
-
-        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,],])
-
-        disp = numpy.array([[ 8.1, 9.1,],
-                            [ 8.2, 9.2,],
-                            [ 8.3, 9.3,],
-                            [ 8.4, 9.4,],
-                            [ 8.2, 9.2,],
-                            [ 8.3, 9.3,],
-                            [-8.6, 9.6,],
-                            [-8.8, 9.8,],])
-
-        if testCase == "slip":
-            dispIncr = numpy.array([[ 9.1, 7.1,],
-                                    [ 9.2, 7.2,],
-                                    [ 9.3, 7.3,],
-                                    [ 9.4, 7.4,],
-                                    [ 9.2, 7.2,],
-                                    [ 9.3, 7.3,],
-                                    [-1.6, 2.6,],
-                                    [-1.8, 2.8,],])            
+        fieldT = numpy.array([[ 8.1, 9.1,],
+                              [ 8.2, 9.2,],
+                              [ 8.3, 9.3,],
+                              [ 8.4, 9.4,],
+                              [ 8.2, 9.2,],
+                              [ 8.3, 9.3,],
+                              [-8.6, 9.6,],
+                              [-8.8, 9.8,],])
+
+        if testCase == "stick":
+            fieldTIncr = numpy.array([[ 1.1, 2.1,],
+                                      [ 1.2, 2.2,],
+                                      [ 1.3, 2.3,],
+                                      [ 1.4, 2.4,],
+                                      [ 1.2, 2.2,],
+                                      [ 1.3, 2.3,],
+                                      [-21.6, 2.6,],
+                                      [-21.8, 2.8,],])            
+        elif testCase == "slip":
+            fieldTIncr = numpy.array([[ 9.1, 7.1,],
+                                      [ 9.2, 7.2,],
+                                      [ 9.3, 7.3,],
+                                      [ 9.4, 7.4,],
+                                      [ 9.2, 7.5,],
+                                      [ 9.3, 7.6,],
+                                      [-1.6, 2.6,],
+                                      [-1.8, 2.8,],])            
         elif testCase == "open":
-            dispIncr = numpy.array([[ 9.1, 7.1,],
-                                    [ 9.2, 7.2,],
-                                    [ 9.3, 7.3,],
-                                    [ 9.4, 7.4,],
-                                    [ 9.2, 7.2,],
-                                    [ 9.3, 7.3,],
-                                    [ +10.6, -10.6,],
-                                    [ +10.8, -10.8,],])            
+            fieldTIncr = numpy.array([[ 9.1, 7.1,],
+                                      [ 9.2, 7.2,],
+                                      [ 9.3, 7.3,],
+                                      [ 9.4, 7.4,],
+                                      [ 9.5, 7.5,],
+                                      [ 9.6, 7.6,],
+                                      [ +10.6, -10.6,],
+                                      [ +10.8, -10.8,],])            
 
 
     elif cell == "tri3d":
@@ -139,22 +131,6 @@ if cell == "tri3" or cell == "tri3d" or cell == "quad4":
                          [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],
-             [-1.0, +6.1, -0.9, -0.8, -0.7, -0.6],
-             [-1.1, -0.9, +6.2, -2.1,  0.0,  0.0],
-             [-1.2, -0.8, -2.1, +6.3,  0.0,  0.0],
-             [-1.3, -0.7,  0.0,  0.0, +6.4, -1.1],
-             [-1.4, -0.6,  0.0,  0.0, -1.1, +6.5]])
-
-        jacobianP = numpy.array(
-            [[+5.0, -1.0, -1.1, -1.2, -1.3, -1.4],
-             [-1.0, +5.1, -0.9, -0.8, -0.7, -0.6],
-             [-1.1, -0.9, +5.2, -2.1,  0.0,  0.0],
-             [-1.2, -0.8, -2.1, +5.3,  0.0,  0.0],
-             [-1.3, -0.7,  0.0,  0.0, +5.4, -1.1],
-             [-1.4, -0.6,  0.0,  0.0, -1.1, +5.5]])
-
         disp = numpy.array([[ 6.1, 8.1,],
                             [ 6.2, 8.2,],
                             [ 6.3, 8.3,],
@@ -218,17 +194,6 @@ if cell == "tri3" or cell == "tri3d" or cell == "quad4":
                          [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,],])
-
-        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,],])
 
         disp = numpy.array([[ 8.1, 9.1,],
                             [ 8.3, 9.3,],
@@ -268,71 +233,31 @@ if cell == "tri3" or cell == "tri3d" or cell == "quad4":
     # ------------------------------------------------------------------
     lagrangeScale = lengthScale**1
 
-    fieldTpdt = fieldT + fieldIncr
-
-    fieldTpdt = globalToFault(fieldTpdt, C)
+    fieldTpdt = (fieldT + fieldTIncr) / lengthScale
+    fieldTpdtFault = globalToFault(fieldTpdt[indexL,:], C)
 
-    tractionShear = abs(fieldTpdt[:,0])
-    tractionNormal = fieldTpdt[:,1]
+    tractionShear = abs(fieldTpdtFault[:,0])
+    tractionNormal = fieldTpdtFault[:,1]
 
     print "tractionShear",tractionShear
     print "tractionNormal",tractionNormal
 
-    friction = -0.6 * tractionNormal;
-
-    print "friction",friction
-
-    dlagrange0 = (friction - tractionShear) * fieldTpdt[:,0] / tractionShear
-  
-    print "dlagrange0",dlagrange0
-
-    if testCase == "slip": 
-        dLagrange = numpy.vstack((dlagrange0, dlagrange1))
-        dLagrange = numpy.transpose(dLagrange)
-        dLagrange = faultToGlobal(dLagrange, C).reshape(m)
-    elif testCase == "open":
-        dLagrange = numpy.reshape(disp+dispIncr, n)
-        dLagrange = -dLagrange[indexL]
+    tractionRheologyFault = numpy.zeros((m, DOF))
+    tractionRheologyFault[:,0] = -0.6 * tractionNormal;
+    tractionRheologyFault[:,1] = tractionNormal
+    print "tractionRheology",tractionRheologyFault
 
-    print "dLagrange \n", dLagrange
+    tractionRheology = faultToGlobal(tractionRheologyFault, C)
+    tractionInternal = fieldTpdt[indexL,:]
+    tractionResidual = tractionRheology - tractionInternal
+    print "tractionResidual",tractionResidual
 
-    L /= lengthScale**1
-    RHS = numpy.dot(numpy.transpose(L),dLagrange)
-    print "RHS",RHS
-    duN = numpy.dot(inv(jacobianN),RHS)
-    duP = -numpy.dot(inv(jacobianP),RHS)
-    
-    dispRelIncr = duP - duN
-
-    dispTpdt = disp + dispIncr
-    dispTpdt = numpy.reshape(dispTpdt, n)
+    residual = numpy.zeros(fieldT.shape)
+    residual[indexN,:] = +tractionResidual
+    residual[indexP,:] = -tractionResidual
+    residual[indexL,:] = tractionInternal * (fieldTpdt[indexP,:] - fieldTpdt[indexN,:]) * lagrangeScale
 
-    slipVertex = dispRelIncr + dispTpdt[indexP]-dispTpdt[indexN]
-    slipVertex = numpy.reshape(slipVertex, (m/DOF,DOF))
-    slipVertex = globalToFault(slipVertex, C)
-    mask = slipVertex[:,1] < 0.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
-
-    dispIncrE = dispIncr
-    dispIncrE = numpy.reshape(dispIncrE, n)
-    dispIncrE[indexL] = dispIncrE[indexL] + dLagrange
-    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))
-    slipVertex = globalToFault(slipVertex, C)
-
-    print "dispIncrE\n", printdata(dispIncrE)
-    print "slipVertexE\n", printdata(slipVertex)
+    print "residual \n",printdata(residual)
 
 
 # ----------------------------------------------------------------------



More information about the CIG-COMMITS mailing list