[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
Wed Nov 5 15:40:38 PST 2014
Repository : https://github.com/geodynamics/pylith
On branch : baagaard/dynrup-new-lagrange
Link : https://github.com/geodynamics/pylith/compare/f33c75b19fd60eedb2a3405db76a1fee333bb1d7...5b6d812b1612809fea3bd331c4e5af98c25a536a
>---------------------------------------------------------------
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