[cig-commits] r7395 - in short/3D/PyLith/trunk: libsrc/faults
unittests/libtests/faults unittests/libtests/faults/data
brad at geodynamics.org
brad at geodynamics.org
Fri Jun 22 15:29:22 PDT 2007
Author: brad
Date: 2007-06-22 15:29:22 -0700 (Fri, 22 Jun 2007)
New Revision: 7395
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc
Log:
Fixed bug in conditioning Jacobian matrix. Only scale the constraint matrix in the upper right half of the Jacobian matrix. Fixed bug in flipping fault orientation- only flip up-dip to correct sense of slip. Updated unit tests accordingly.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-06-22 22:23:40 UTC (rev 7394)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-06-22 22:29:22 UTC (rev 7395)
@@ -233,8 +233,7 @@
if (2 == cohesiveDim) {
// Check orientation of first vertex, if dot product of fault
- // normal with preferred normal is negative, flip all along-strike
- // and fault-normal directions (keep up-dip direction the same).
+ // normal with preferred normal is negative, flip up/down dip direction.
// If the user gives the correct normal direction, we should end
// up with left-lateral-slip, reverse-slip, and fault-opening for
// positive slip values.
@@ -254,15 +253,15 @@
const real_section_type::value_type* vertexOrient =
_orientation->restrictPoint(*v_iter);
assert(9 == _orientation->getFiberDimension(*v_iter));
- // Flip along-strike direction
+ // Keep along-strike direction
for (int iDim=0; iDim < 3; ++iDim)
+ vertexDir[iDim] = vertexOrient[iDim];
+ // Flip up-dip direction
+ for (int iDim=3; iDim < 6; ++iDim)
vertexDir[iDim] = -vertexOrient[iDim];
- // Keep up-dip direction
- for (int iDim=3; iDim < 6; ++iDim)
+ // Keep normal direction
+ for (int iDim=6; iDim < 9; ++iDim)
vertexDir[iDim] = vertexOrient[iDim];
- // Flip normal direction
- for (int iDim=6; iDim < 9; ++iDim)
- vertexDir[iDim] = -vertexOrient[iDim];
// Update direction
_orientation->updatePoint(*v_iter, &vertexDir[0]);
@@ -350,6 +349,7 @@
const double density = matprops[0];
const double vs = matprops[1];
const double mu = density * vs*vs;
+ //const double mu = 1.0;
_pseudoStiffness->updatePoint(*v_iter, &mu);
} // for
} // initialize
@@ -452,7 +452,7 @@
// constraint vertices (k) of the cohesive cells
for (int iDim=0; iDim < spaceDim; ++iDim)
cellResidual[indexK*spaceDim+iDim] =
- cellSlip[iConstraint*spaceDim+iDim] * pseudoStiffness;
+ cellSlip[iConstraint*spaceDim+iDim];
// Get orientation at constraint vertex
const real_section_type::value_type* constraintOrient =
@@ -570,7 +570,7 @@
cellMatrix[row*numCorners*spaceDim+col] =
-constraintOrient[kDim*spaceDim+iDim]*pseudoStiffness;
cellMatrix[col*numCorners*spaceDim+row] =
- cellMatrix[row*numCorners*spaceDim+col]; // symmetric
+ -constraintOrient[kDim*spaceDim+iDim];
} // for
// Entries associated with constraint forces applied at node j
@@ -581,7 +581,7 @@
cellMatrix[row*numCorners*spaceDim+col] =
constraintOrient[kDim*spaceDim+jDim]*pseudoStiffness;
cellMatrix[col*numCorners*spaceDim+row] =
- cellMatrix[row*numCorners*spaceDim+col]; // symmetric
+ constraintOrient[kDim*spaceDim+jDim];
} // for
} // for
err = PetscLogFlops(numConstraintVert*spaceDim*spaceDim*4);
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-06-22 22:23:40 UTC (rev 7394)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-06-22 22:29:22 UTC (rev 7395)
@@ -23,6 +23,15 @@
* The ordering of vertices in a cohesive cell is the vertices on the
* one side of the fault, the corresponding entries on the other side
* of the fault, and then the corresponding constraint vertices.
+ *
+ * [ K aC^T ] [ U ] = [ Fe ]
+ * [ C 0 ] [ Fi ] = [ D ]
+ *
+ * where K is the stiffness matrix, C is the matrix of Lagrange
+ * constraints, U is the displacement field, Fe is the vector of
+ * external forces, Fi is the vector of Lagrange multipers (forces), D
+ * is the fault slip, and "a" is the conditioning value (taken to be
+ * the shear modulus).
*/
#if !defined(pylith_faults_faultcohesivekin_hh)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2007-06-22 22:23:40 UTC (rev 7394)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2007-06-22 22:29:22 UTC (rev 7395)
@@ -235,9 +235,14 @@
CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
const real_section_type::value_type* vals =
residual->restrictPoint(*v_iter);
+
+ const bool isConstraint =
+ (fault._constraintVert.end() != fault._constraintVert.find(*v_iter));
+ const double pseudoStiffness =
+ (!isConstraint) ? _data->pseudoStiffness : 1.0;
for (int i=0; i < fiberDimE; ++i) {
const int index = iVertex*spaceDim+i;
- const double valE = valsE[index] * _data->pseudoStiffness;
+ const double valE = valsE[index] * pseudoStiffness;
if (valE > tolerance)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
else
@@ -325,7 +330,9 @@
for (int iRow=0; iRow < nrows; ++iRow)
for (int iCol=0; iCol < ncols; ++iCol) {
const int index = ncols*iRow+iCol;
- const double valE = valsE[index] * _data->pseudoStiffness;
+ const double pseudoStiffness =
+ (iRow > iCol) ? 1.0 : _data->pseudoStiffness;
+ const double valE = valsE[index] * pseudoStiffness;
#if 0 // DEBUGGING
if (fabs(valE-vals[index]) > tolerance)
std::cout << "ERROR: iRow: " << iRow << ", iCol: " << iCol
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc 2007-06-22 22:23:40 UTC (rev 7394)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc 2007-06-22 22:29:22 UTC (rev 7395)
@@ -97,9 +97,9 @@
const int pylith::faults::CohesiveKinDataTet4f::_numConstraintVert = 3;
const double pylith::faults::CohesiveKinDataTet4f::_orientation[] = {
- 0.0, +1.0, 0.0, 0.0, 0.0, +1.0, +1.0, 0.0, 0.0,
- 0.0, +1.0, 0.0, 0.0, 0.0, +1.0, +1.0, 0.0, 0.0,
- 0.0, +1.0, 0.0, 0.0, 0.0, +1.0, +1.0, 0.0, 0.0,
+ 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, -1.0, 0.0, 0.0,
+ 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, -1.0, 0.0, 0.0,
+ 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, -1.0, 0.0, 0.0,
};
const int pylith::faults::CohesiveKinDataTet4f::_constraintVertices[] = {
@@ -112,15 +112,15 @@
const double pylith::faults::CohesiveKinDataTet4f::_valsResidual[] = {
0.0, 0.0, 0.0,
- 9.7, 7.7, 8.7, // 3
- 9.9, 7.9, 8.9, // 4
- 9.1, 7.1, 8.1, // 5
+ -9.7, -7.7, -8.7, // 3
+ -9.9, -7.9, -8.9, // 4
+ -9.1, -7.1, -8.1, // 5
0.0, 0.0, 0.0,
- -9.7, -7.7, -8.7, // 7
+ +9.7, +7.7, +8.7, // 7
1.07974939836, -0.32861938211, 0.04694562602, // 8
- -9.9, -7.9, -8.9, // 4
+ +9.9, +7.9, +8.9, // 4
1.00381374723, -0.33460458241, 0.08365114560, // 10
- -9.1, -7.1, -8.1, // 5
+ +9.1, +7.1, +8.1, // 5
0.90493237602, -0.32577565537, 0.10859188512, // 12
};
@@ -164,7 +164,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,-1.0, // 8
+ 0.0, 0.0,+1.0, // 8
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -175,7 +175,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- -1.0, 0.0, 0.0, // 8
+ +1.0, 0.0, 0.0, // 8
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -186,7 +186,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,-1.0, 0.0, // 8
+ 0.0,+1.0, 0.0, // 8
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -199,7 +199,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,-1.0, // 10
+ 0.0, 0.0,+1.0, // 10
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 4y
@@ -210,7 +210,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- -1.0, 0.0, 0.0, // 10
+ +1.0, 0.0, 0.0, // 10
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 4z
@@ -221,7 +221,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,-1.0, 0.0, // 10
+ 0.0,+1.0, 0.0, // 10
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 5x
@@ -234,7 +234,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,-1.0, // 12
+ 0.0, 0.0,+1.0, // 12
0.0, 0.0, 0.0, // 5y
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -245,7 +245,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- -1.0, 0.0, 0.0, // 12
+ +1.0, 0.0, 0.0, // 12
0.0, 0.0, 0.0, // 5z
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -256,7 +256,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,-1.0, 0.0, // 12
+ 0.0,+1.0, 0.0, // 12
0.0, 0.0, 0.0, // 6x
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -296,7 +296,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,+1.0, // 8
+ 0.0, 0.0,-1.0, // 8
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -307,7 +307,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- +1.0, 0.0, 0.0, // 8
+ -1.0, 0.0, 0.0, // 8
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -318,39 +318,39 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,+1.0, 0.0, // 8
+ 0.0,-1.0, 0.0, // 8
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 8x
- 0.0,-1.0, 0.0, // 3
+ 0.0,+1.0, 0.0, // 3
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,+1.0, 0.0, // 7
+ 0.0,-1.0, 0.0, // 7
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 8y
- 0.0, 0.0,-1.0, // 3
+ 0.0, 0.0,+1.0, // 3
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,+1.0, // 7
+ 0.0, 0.0,-1.0, // 7
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 8z
- -1.0, 0.0, 0.0, // 3
+ +1.0, 0.0, 0.0, // 3
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- +1.0, 0.0, 0.0, // 7
+ -1.0, 0.0, 0.0, // 7
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -364,7 +364,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,+1.0, // 10
+ 0.0, 0.0,-1.0, // 10
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 9y
@@ -375,7 +375,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- +1.0, 0.0, 0.0, // 10
+ -1.0, 0.0, 0.0, // 10
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 9z
@@ -386,39 +386,39 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,+1.0, 0.0, // 10
+ 0.0,-1.0, 0.0, // 10
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 10x
0.0, 0.0, 0.0,
- 0.0,-1.0, 0.0, // 4
+ 0.0,+1.0, 0.0, // 4
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,+1.0, 0.0, // 9
+ 0.0,-1.0, 0.0, // 9
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 10y
0.0, 0.0, 0.0,
- 0.0, 0.0,-1.0, // 4
+ 0.0, 0.0,+1.0, // 4
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,+1.0, // 9
+ 0.0, 0.0,-1.0, // 9
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 10z
0.0, 0.0, 0.0,
- -1.0, 0.0, 0.0, // 4
+ +1.0, 0.0, 0.0, // 4
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- +1.0, 0.0, 0.0, // 9
+ -1.0, 0.0, 0.0, // 9
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -432,7 +432,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,+1.0, // 12
+ 0.0, 0.0,-1.0, // 12
0.0, 0.0, 0.0, // 11y
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -443,7 +443,7 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- +1.0, 0.0, 0.0, // 12
+ -1.0, 0.0, 0.0, // 12
0.0, 0.0, 0.0, // 11z
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
@@ -454,39 +454,39 @@
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,+1.0, 0.0, // 12
+ 0.0,-1.0, 0.0, // 12
0.0, 0.0, 0.0, // 12x
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,-1.0, 0.0, // 5
+ 0.0,+1.0, 0.0, // 5
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0,+1.0, 0.0, // 11
+ 0.0,-1.0, 0.0, // 11
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 12y
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,-1.0, // 5
+ 0.0, 0.0,+1.0, // 5
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- 0.0, 0.0,+1.0, // 11
+ 0.0, 0.0,-1.0, // 11
0.0, 0.0, 0.0,
0.0, 0.0, 0.0, // 12z
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- -1.0, 0.0, 0.0, // 5
+ +1.0, 0.0, 0.0, // 5
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
- +1.0, 0.0, 0.0, // 11
+ -1.0, 0.0, 0.0, // 11
0.0, 0.0, 0.0,
};
More information about the cig-commits
mailing list