[cig-commits] r18953 - in short/3D/PyLith/branches/v1.6-revisedfault: libsrc/pylith/faults unittests/libtests/faults unittests/libtests/faults/data
brad at geodynamics.org
brad at geodynamics.org
Tue Sep 20 18:00:49 PDT 2011
Author: brad
Date: 2011-09-20 18:00:48 -0700 (Tue, 20 Sep 2011)
New Revision: 18953
Modified:
short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
Log:
Started debugging revised fault implementation (fixed some bugs).
Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-09-20 19:50:22 UTC (rev 18952)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-09-21 01:00:48 UTC (rev 18953)
@@ -344,7 +344,8 @@
const int iB = iBasis*spaceDim;
for (int iDim = 0; iDim < spaceDim; ++iDim)
for (int kDim = 0; kDim < spaceDim; ++kDim)
- slipGlobalCell[iB+iDim] += slipCell[iB+kDim] * orientationCell[iB+kDim];
+ slipGlobalCell[iB+iDim] += slipCell[iB+kDim] *
+ orientationCell[iB*spaceDim + kDim*spaceDim + iDim];
} // for
// Compute action for positive side of fault and Lagrange constraints
@@ -353,32 +354,49 @@
const int iQ = iQuad * numBasis;
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQ+iBasis];
+
+ // Add entries to residual at
+ // iP = DOF at positive node
+ // iN = DOF at negative node
+ // iL = DOF at constraint node
+
+ // Indices for positive vertex
+ const int iBP = 0*numBasis*spaceDim + iBasis*spaceDim;
+
+ // Indices for negative vertex
+ const int iBN = 1*numBasis*spaceDim + iBasis*spaceDim;
+
+ // Indices for Lagrange vertex
+ const int iBL = 2*numBasis*spaceDim + iBasis*spaceDim;
+
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
const double valIJ = valI * basis[iQ+jBasis];
- // positive side of the fault
+ // Indices for positive vertex
+ const int jBP = 0*numBasis*spaceDim + jBasis*spaceDim;
+
+ // Indices for negative vertex
+ const int jBN = 1*numBasis*spaceDim + jBasis*spaceDim;
+
+ // Indices for Lagrange vertex
+ const int jBL = 2*numBasis*spaceDim + jBasis*spaceDim;
+
for (int iDim=0; iDim < spaceDim; ++iDim) {
- residualCell[0*numBasis*spaceDim+iBasis*spaceDim+iDim] +=
- dispTpdtCell[2*numBasis*spaceDim+jBasis*spaceDim+iDim] * valIJ;
+ // positive side of the fault
+ residualCell[iBP + iDim] += dispTpdtCell[jBL + iDim] * valIJ;
- // Lagrange constraints
- residualCell[2*numBasis*spaceDim+iBasis*spaceDim+iDim] +=
- slipGlobalCell[jBasis*spaceDim+iDim] -
- (dispTpdtCell[1*numBasis*spaceDim+jBasis*spaceDim+iDim] -
- dispTpdtCell[0*numBasis*spaceDim+jBasis*spaceDim+iDim])
- * valIJ;
- } // forg
+ // negative side of the fault
+ residualCell[iBN + iDim] = -residualCell[iBP + iDim];
+
+ // Lagrange constraint
+ residualCell[iBL + iDim] +=
+ valIJ * (slipGlobalCell[jBasis*spaceDim+iDim]
+ - dispTpdtCell[jBP + iDim] + dispTpdtCell[jBN + iDim]);
+ } // for
} // for
} // for
} // for
- // Compute action for negative side of fault
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- for (int iDim=0; iDim < spaceDim; ++iDim)
- residualCell[1*numBasis*spaceDim+iBasis*spaceDim+iDim] =
- -residualCell[0*numBasis*spaceDim+iBasis*spaceDim+iDim];
- } // for
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
@@ -510,33 +528,47 @@
jacobianCell = 0.0;
- // Compute Jacobian for constrants, positive side
+ const int rowSize = 3*numBasis*spaceDim;
+
+ // Compute Jacobian for constrants
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQ+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
const double valIJ = valI * basis[iQ+jBasis];
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- const int iBlock = (0*numBasis*spaceDim + iBasis*spaceDim + iDim) *
- (3*numBasis*spaceDim);
- const int jBlock = (2*numBasis*spaceDim + jBasis*spaceDim + iDim);
- jacobianCell[iBlock+jBlock] += valIJ;
+
+ // First index for positive, negative, and Lagrange vertices
+ const int iBP = 0*numBasis*spaceDim + iBasis*spaceDim;
+ const int iBN = 1*numBasis*spaceDim + iBasis*spaceDim;
+ const int iBL = 2*numBasis*spaceDim + iBasis*spaceDim;
+
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ // Add entries to Jacobian at (i,j) where
+ // iP = DOF at positive node, jL = DOF at constraint node
+ // iL = DOF at constraint node, jP = DOF at positive node
+ // iN = DOF at negative node, jL = DOF at constraint node
+ // iL = DOF at constraint node, jN = DOF at negative node
+
+ const int iP = (iBP + iDim) * rowSize; // row
+ const int jP = (iBP + iDim); // col
+
+ // Indices for negative vertex
+ const int iN = (iBN + iDim) * rowSize; // row
+ const int jN = (iBN + iDim); // col
+
+ // Indices for Lagrange vertex
+ const int iL = (iBL + iDim) * rowSize; // row
+ const int jL = (iBL + iDim); // col
+
+ jacobianCell[iP + jL] -= valIJ;
+ jacobianCell[iL + jP] -= valIJ;
+
+ jacobianCell[iN + jL] += valIJ;
+ jacobianCell[iL + jN] += valIJ;
} // for
} // for
} // for
-
- // Compute Jacobian for constrants, negative side
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- const int iBlock = (0*numBasis*spaceDim + iBasis*spaceDim + iDim) *
- (3*numBasis*spaceDim);
- const int jBlock = (2*numBasis*spaceDim + jBasis*spaceDim + iDim);
- jacobianCell[iBlock+jBlock] = -jacobianCell[iBlock+jBlock];
- } // for
- } // for
- } // for
} // for
#if defined(DETAILED_EVENT_LOGGING)
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveKin.cc 2011-09-20 19:50:22 UTC (rev 18952)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveKin.cc 2011-09-21 01:00:48 UTC (rev 18953)
@@ -264,7 +264,7 @@
fault.useSolnIncr(false);
fault.integrateResidual(residual, t, &fields);
- //residual.view("RESIDUAL"); // DEBUGGING
+ residual.view("RESIDUAL"); // DEBUGGING
// Check values
const double* valsE = _data->residual;
@@ -295,7 +295,7 @@
fault.useSolnIncr(true);
fault.integrateResidual(residual, t, &fields);
- //residual.view("RESIDUAL"); // DEBUGGING
+ residual.view("RESIDUAL"); // DEBUGGING
// Check values
const double* valsE = _data->residualIncr;
@@ -356,7 +356,7 @@
jacobian.assemble("final_assembly");
- //MatView(jacobian.matrix(), PETSC_VIEWER_STDOUT_WORLD); // DEBUGGING
+ MatView(jacobian.matrix(), PETSC_VIEWER_STDOUT_WORLD); // DEBUGGING
const double* valsE = _data->jacobian;
const int nrowsE = dispSection->sizeWithBC();
@@ -388,7 +388,7 @@
for (int iCol=0; iCol < ncols; ++iCol) {
const int index = ncols*iRow+iCol;
const double valE = valsE[index];
-#if 0 // DEBUGGING
+#if 1 // DEBUGGING
if (fabs(valE-vals[index]) > tolerance)
std::cout << "ERROR: iRow: " << iRow << ", iCol: " << iCol
<< "valE: " << valE
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2011-09-20 19:50:22 UTC (rev 18952)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2011-09-21 01:00:48 UTC (rev 18953)
@@ -133,7 +133,7 @@
7.5,
0.0,
-7.5,
- -0.2+1.89546413727,
+ -0.2+1.89546413727,
};
const double pylith::faults::CohesiveKinDataLine2::_residual[] = {
@@ -141,7 +141,7 @@
7.5, // 3
0.0,
-7.5, // 5
- -0.2+1.89546413727,
+ -0.2+1.89546413727,
};
const double pylith::faults::CohesiveKinDataLine2::_jacobian[] = {
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2011-09-20 19:50:22 UTC (rev 18952)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2011-09-21 01:00:48 UTC (rev 18953)
@@ -166,24 +166,26 @@
const double pylith::faults::CohesiveKinDataTri3::_residual[] = {
0.0, 0.0,
- -9.6, -8.6, // 3
- -9.8, -8.8, // 4
+ -8.7, -9.7, // 3
+ -8.7, -9.7, // 4
0.0, 0.0,
- +9.6, +8.6, // 6
- +9.8, +8.8, // 7
- 0.3+1.89546413727, 0.3+0.08241148423, // 8
- 0.4+1.77538035254, 0.4+0.14794836271, // 9
+ +8.7, +9.7, // 6
+ +8.7, +9.7, // 7
+ 0.5*(0.08241148423+0.14794836271) - 0.5*(3.5-3.2 + 3.7-3.3),
+ 0.5*(1.89546413727+1.77538035254) - 0.5*(4.5-4.2 + 4.7-4.3), // 8
+ 0.5*(0.08241148423+0.14794836271) - 0.5*(3.5-3.2 + 3.7-3.3),
+ 0.5*(1.89546413727+1.77538035254) - 0.5*(4.5-4.2 + 4.7-4.3), // 9
};
const double pylith::faults::CohesiveKinDataTri3::_residualIncr[] = {
0.0, 0.0,
- -9.6, -8.6, // 3
- -9.8, -8.8, // 4
+ -8.6, -9.6, // 3
+ -8.8, -9.8, // 4
0.0, 0.0,
- +9.6, +8.6, // 6
- +9.8, +8.8, // 7
- 0.3+1.89546413727, 0.3+0.08241148423, // 8
- 0.4+1.77538035254, 0.4+0.14794836271, // 9
+ +8.6, +9.6, // 6
+ +8.8, +9.8, // 7
+ 0.3+0.08241148423, 0.3+1.89546413727, // 8
+ 0.4+0.14794836271, 0.4+1.77538035254, // 9
};
const double pylith::faults::CohesiveKinDataTri3::_jacobian[] = {
More information about the CIG-COMMITS
mailing list