[cig-commits] r16308 - in short/3D/PyLith/trunk: libsrc/faults unittests/libtests/faults unittests/libtests/faults/data

brad at geodynamics.org brad at geodynamics.org
Mon Feb 22 13:03:01 PST 2010


Author: brad
Date: 2010-02-22 13:03:01 -0800 (Mon, 22 Feb 2010)
New Revision: 16308

Added:
   short/3D/PyLith/trunk/unittests/libtests/faults/data/cohesivedyn.py
Modified:
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
Log:
Completed updating unit tests for weighted average of initial tractions. Removed some debugging output.

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc	2010-02-22 19:44:40 UTC (rev 16307)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc	2010-02-22 21:03:01 UTC (rev 16308)
@@ -1088,7 +1088,7 @@
     tractionSection->updatePoint(*v_iter, &tractionVertex[0]);
   } // for
 
-  traction.view("INITIAL TRACTIONS"); // DEBUGGING
+  //traction.view("INITIAL TRACTIONS"); // DEBUGGING
 } // _getInitialTractions
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc	2010-02-22 19:44:40 UTC (rev 16307)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc	2010-02-22 21:03:01 UTC (rev 16308)
@@ -145,6 +145,7 @@
 
   // Check orientation
   //fault._fields->get("orientation").view("ORIENTATION"); // DEBUGGING
+
   const ALE::Obj<RealSection>& orientationSection = 
     fault._fields->get("orientation").section();
   CPPUNIT_ASSERT(!orientationSection.isNull());
@@ -375,9 +376,11 @@
       for (int i = 0; i < fiberDimE; ++i) {
         const int index = iVertex * spaceDim + i;
         const double valE = valsE[index];
+#if 0 // DEBUGGING
 	std::cout << "valE: " << valE
 		  << ", val: " << vals[i]
 		  << std::endl;
+#endif // DEBUGGING
         if (fabs(valE) > tolerance)
           CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
         else
@@ -424,9 +427,11 @@
       for (int i = 0; i < fiberDimE; ++i) {
         const int index = iVertex * spaceDim + i;
         const double valE = valsE[index];
+#if 0 // DEBUGGING
 	std::cout << "valE: " << valE
 		  << ", val: " << vals[i]
 		  << std::endl;
+#endif // DEBUGGING
         if (fabs(valE) > tolerance)
           CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
         else
@@ -494,9 +499,11 @@
       for (int i = 0; i < fiberDimE; ++i) {
         const int index = iVertex * spaceDim + i;
         const double valE = valsE[index];
+#if 0 // DEBUGGING
 	std::cout << "valE: " << valE
 		  << ", val: " << vals[i]
 		  << std::endl;
+#endif // DEBUGGING
         if (fabs(valE) > tolerance)
           CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
         else
@@ -543,9 +550,11 @@
       for (int i = 0; i < fiberDimE; ++i) {
         const int index = iVertex * spaceDim + i;
         const double valE = valsE[index];
+#if 0 // DEBUGGING
 	std::cout << "valE: " << valE
 		  << ", val: " << vals[i]
 		  << std::endl;
+#endif // DEBUGGING
         if (fabs(valE) > tolerance)
           CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
         else

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2010-02-22 19:44:40 UTC (rev 16307)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2010-02-22 21:03:01 UTC (rev 16308)
@@ -332,7 +332,7 @@
     fault.useSolnIncr(false);
     fault.integrateResidualAssembled(residual, t, &fields);
 
-    residual.view("RESIDUAL"); // DEBUGGING
+    //residual.view("RESIDUAL"); // DEBUGGING
 
     // Check values
     const double* valsE = _data->residual;
@@ -574,9 +574,7 @@
   CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
   jacobian.complete();
 
-#if 0 // DEBUGGING
-  jacobian.view("JACOBIAN");
-#endif // DEBUGGING
+  // jacobian.view("JACOBIAN"); // DEBUGGING
 
   const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
   CPPUNIT_ASSERT(!jacobianSection.isNull());

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc	2010-02-22 19:44:40 UTC (rev 16307)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc	2010-02-22 21:03:01 UTC (rev 16308)
@@ -1427,18 +1427,18 @@
   6.5, 8.5, 10.5, // 15
   6.7, 8.7, 10.7, // 16
   6.9, 8.9, 10.9, // 17
-  -3.963897477671262,  -5.477083063322538,  -10.4, // 18
-  -4.120121221785221,  -5.634915064993491,  -10.6, // 19
-  -4.276416812742042,  -5.792739610437928,  -10.8, // 20 
-  -3.474281747841006,  -4.919475110692280,  -10.0, // 21
+  -3.94113616, -5.44660605,  -10.4, // 18
+  -4.11249353, -5.62478649,  -10.6, // 19
+  -4.28408279, -5.80283922,  -10.8, // 20 
+  -3.49677632, -4.95014952,  -10.0, // 21
 };
 
 // :TODO: Update slip values based on changes in Lagrange multiplier values
 const double pylith::faults::CohesiveDynDataHex8::_slipSlipE[] = {
-  20.727794955342528,  27.754166126645078,                   0.0,
-  21.440242443570444,  28.469830129986981,                   0.0,
-  22.152833625484085,  29.185479220875859,                   0.0,
-  18.948563495682009,  25.838950221384557,                   0.0,
+  20.68227232,  27.69321209,  0.0,
+  21.42498706,  28.44957298,  0.0,
+  22.16816557,  29.20567845,  0.0,
+  18.99355263,  25.90029904,  0.0,
 };
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc	2010-02-22 19:44:40 UTC (rev 16307)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc	2010-02-22 21:03:01 UTC (rev 16308)
@@ -532,16 +532,16 @@
   8.6, 9.6, 10.6, // 7
   8.8, 9.8, 10.8, // 8
   8.0, 9.0, 10.0, // 9
-  -6.90156, -7.80418, -10.7, // 10
-  -7.08708, -7.9903, -10.9, // 11
-  -6.27887, -7.17083, -10.1, // 12
+  -6.88824796, -7.78925381, -10.7, // 10
+  -7.0870762, -7.99029956, -10.9, // 11
+  -6.29211869, -7.18581852, -10.1, // 12
 };
 
 // :TODO: Update slip values based on changes in Lagrange multiplier values
 const double pylith::faults::CohesiveDynDataTet4::_slipSlipE[] = {
-  31.203110740589327,  35.008368147978267,                   0.0,
-  31.974152400186210,  35.780599114494088,                   0.0,
-  28.757749464750159,  32.541663868006758,                   0.0,
+  31.17649592,  34.97850762,  0.0,
+  31.9741524,  35.78059911,  0.0,
+  28.78423738,  32.57163703,  0.0,
 };
 
 // ----------------------------------------------------------------------

Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/cohesivedyn.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/cohesivedyn.py	                        (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/cohesivedyn.py	2010-02-22 21:03:01 UTC (rev 16308)
@@ -0,0 +1,103 @@
+cell = "hex8"
+dim = "3d"
+
+import numpy
+
+# ----------------------------------------------------------------------
+if dim == "3d":
+    if cell == "tet4":
+
+        fieldT = numpy.array([[7.7, 8.7, 9.7],
+                              [7.9, 8.9, 9.9],
+                              [7.1, 8.1, 9.1]])
+        fieldIncr = numpy.array([[8.7, 9.7, -10.7],
+                                 [8.9, 9.9, -10.9],
+                                 [8.1, 9.1, -10.1]])
+        initialTract = numpy.array([[1.1, 2.1, -3.1],
+                                    [1.1, 2.1, -3.1],
+                                    [1.1, 2.1, -3.1]])
+        area = numpy.array([1.0/3.0, 1.0/3.0, 1.0/3.0])
+    
+
+    elif cell == "hex8":
+        fieldT = numpy.array([[5.4, 7.4, 9.4],
+                              [5.6, 7.6, 9.6],
+                              [5.8, 7.8, 9.8],
+                              [5.0, 7.0, 9.0]])
+        fieldIncr = numpy.array([[6.4, 8.4, -10.4],
+                                 [6.6, 8.6, -10.6],
+                                 [6.8, 8.8, -10.8],
+                                 [6.0, 8.0, -10.0]])
+        initialTract = numpy.array([[1.063397471, 2.063397471, -3.063397471], 
+                                    [1.121132498, 2.121132498, -3.121132498], 
+                                    [1.178867525, 2.178867525, -3.178867525],
+                                    [1.236602552, 2.236602552, -3.236602552]])
+        area = numpy.array([1.0, 1.0, 1.0, 1.0])
+    
+
+        tp = numpy.array([9.4, 9.6, 9.8, 9])
+        tdt = numpy.array([-10.4, -10.6, -10.8, -10])
+        it = numpy.array([-3.0634, -3.12113, -3.17887, -3.2366])
+        
+        tps1 = numpy.array([5.4, 5.6, 5.8, 5])
+        tps2 = numpy.array([7.4, 7.6, 7.8, 7])
+        tdts1 = numpy.array([6.4, 6.6, 6.8, 6])
+        tdts2 = numpy.array([8.4, 8.6, 8.8, 8])
+        its1 = numpy.array([1.0634, 1.12113, 1.17887, 1.2366])
+        its2 = numpy.array([2.0634, 2.12113, 2.17887, 2.2366])
+        
+        tpdts1 = tps1 + tdts1
+        tpdts2 = tps2 + tdts2
+        
+        ts = (tpdts1**2 + tpdts2**2)**0.5
+        nt = tp + tdt + it
+        
+        friction = -0.6 * nt
+        
+        dL0 = (ts - friction) * tpdts1 / ts;
+        dL1 = (ts - friction) * tpdts2 / ts;
+        
+        slipVertex0 = 2 * dL0;
+        slipVertex1 = 2 * dL1;
+
+        tpdts1 = friction * tpdts1 / ts
+        tpdts2 = friction * tpdts2 / ts
+    
+
+    # ------------------------------------------------------------------
+    fieldTpdt = fieldT + fieldIncr
+
+    tractionShear = (fieldTpdt[:,0]**2 + fieldTpdt[:,1]**2)**0.5 / area
+    tractionNormal = initialTract[:,2] + fieldTpdt[:,2] / area
+
+    print "tractionShear",tractionShear
+    print "tractionNormal",tractionNormal
+
+    friction = -0.6 * tractionNormal;
+
+    print "friction",friction
+
+    lagrangeTpdt0 = friction * fieldTpdt[:,0] / tractionShear
+    lagrangeTpdt1 = friction * fieldTpdt[:,1] / tractionShear
+
+    lagrangeIncr0 = lagrangeTpdt0 - fieldT[:,0]
+    lagrangeIncr1 = lagrangeTpdt1 - fieldT[:,1]
+
+    print "lagrangeIncr0",lagrangeIncr0
+    print "lagrangeIncr1",lagrangeIncr1
+
+    dlagrange0 = (tractionShear - friction) * fieldTpdt[:,0] / tractionShear
+    dlagrange1 = (tractionShear - friction) * fieldTpdt[:,1] / tractionShear
+    
+    print "dlagrange0",dlagrange0
+    print "dlagrange1",dlagrange1
+
+    slipVertex0 = 2 * dlagrange0
+    slipVertex1 = 2 * dlagrange1
+
+    print "slipVertex0",slipVertex0
+    print "slipVertex1",slipVertex1
+
+    print "fieldIncr",numpy.array([lagrangeIncr0, lagrangeIncr1]).transpose()
+    print "slip",numpy.array([slipVertex0, slipVertex1]).transpose()
+



More information about the CIG-COMMITS mailing list