[cig-commits] r7122 - in short/3D/PyLith/trunk: libsrc/faults
unittests/libtests/faults unittests/libtests/faults/data
brad at geodynamics.org
brad at geodynamics.org
Sat Jun 9 16:32:52 PDT 2007
Author: brad
Date: 2007-06-09 16:32:51 -0700 (Sat, 09 Jun 2007)
New Revision: 7122
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
Log:
Fixed more bugs in implementing cohesive cells for kinematic earthquake source. Unit tests pass in 1-D. Still need to finish setting test data for some tri3 tests.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-06-09 19:28:28 UTC (rev 7121)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-06-09 23:32:51 UTC (rev 7122)
@@ -89,8 +89,8 @@
const double j2 = jacobian[1];
(*orientation)[0] = j1;
(*orientation)[1] = j2;
- (*orientation)[2] = j2;
- (*orientation)[3] = -j1;
+ (*orientation)[2] = -j2;
+ (*orientation)[3] = j1;
} // _orient2D
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-06-09 19:28:28 UTC (rev 7121)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-06-09 23:32:51 UTC (rev 7122)
@@ -230,8 +230,9 @@
assert(0 != fields);
assert(!mesh.isNull());
- // Subtract constraint forces (which are disp at the constraint
- // DOF) to residual; contributions are at DOF of normal vertices (i and j)
+ // Subtract constraint forces (which are disp at the constraint DOF)
+ // to residual; contributions are at DOF of non-constraint vertices
+ // (i and j) of the cohesive cells
// Get cohesive cells
const ALE::Obj<ALE::Mesh::label_sequence>& cellsCohesive =
@@ -246,25 +247,31 @@
const ALE::Obj<real_section_type>& disp = fields->getHistoryItem(1);
assert(!disp.isNull());
- // Allocate vector for cell values (if necessary)
- _initCellVector();
+ // Allocate vector for cell values
+ const int spaceDim = _quadrature->spaceDim();
+ const int numConstraintVert = _quadrature->numBasis();
+ const int numCorners = 3*numConstraintVert; // cohesive cell
+ double_array cellVector(numCorners*spaceDim);
// Loop over cohesive cells
- const int numConstraintVert = _quadrature->numBasis();
for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
c_iter != cellsCohesiveEnd;
++c_iter) {
- _resetCellVector();
+ cellVector = 0.0;
// Get values at vertices (want constraint forces in disp vector)
const real_section_type::value_type* cellDisp =
mesh->restrict(disp, *c_iter);
// Transfer constraint forces to cell's constribution to residual vector
- for (int i=0; i < numConstraintVert; ++i) {
- const double constraintForce = cellDisp[2*numConstraintVert+i];
- _cellVector[ i] = -constraintForce;
- _cellVector[numConstraintVert+i] = -constraintForce;
+ for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint)
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const int indexI = iConstraint;
+ const int indexJ = iConstraint + numConstraintVert;
+ const int indexK = iConstraint + 2*numConstraintVert;
+ const double constraintForce = cellDisp[indexK*spaceDim+iDim];
+ cellVector[indexI*spaceDim+iDim] = -constraintForce;
+ cellVector[indexJ*spaceDim+iDim] = +constraintForce;
} // for
PetscErrorCode err =
PetscLogFlops(numConstraintVert*2);
@@ -272,7 +279,7 @@
throw std::runtime_error("Logging PETSc flops failed.");
// Update residual (replace, do not add)
- mesh->update(residual, *c_iter, _cellVector);
+ mesh->update(residual, *c_iter, &cellVector[0]);
} // for
} // integrateResidual
@@ -327,9 +334,9 @@
const int jBasis = 3*iConstraint + 1;
const int kBasis = 3*iConstraint + 2;
- // Get orientations at constraint vertex (3rd vertex)
+ // Get orientations at constraint vertex
const real_section_type::value_type* constraintOrient =
- &cellOrientation[kBasis*orientationSize];
+ &cellOrientation[iConstraint*orientationSize];
// Entries associated with constraint forces applied at node i
for (int iDim=0; iDim < spaceDim; ++iDim)
@@ -337,7 +344,7 @@
const int row = iBasis*spaceDim+iDim;
const int col = kBasis*spaceDim+kDim;
cellMatrix[row*numCorners*spaceDim+col] =
- constraintOrient[iDim*spaceDim+kDim];
+ -constraintOrient[iDim*spaceDim+kDim];
cellMatrix[col*numCorners*spaceDim+row] =
cellMatrix[row*numCorners*spaceDim+col]; // symmetric
} // for
@@ -348,7 +355,7 @@
const int row = jBasis*spaceDim+jDim;
const int col = kBasis*spaceDim+kDim;
cellMatrix[row*numCorners*spaceDim+col] =
- -constraintOrient[jDim*spaceDim+kDim];
+ constraintOrient[jDim*spaceDim+kDim];
cellMatrix[col*numCorners*spaceDim+row] =
cellMatrix[row*numCorners*spaceDim+col]; // symmetric
} // for
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2007-06-09 19:28:28 UTC (rev 7121)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2007-06-09 23:32:51 UTC (rev 7122)
@@ -175,8 +175,8 @@
-0.5, 1.0
};
const double orientationE[] = {
- -1.0, 2.0, 2.0, 1.0,
- -0.5, 1.0, 1.0, 0.5
+ -1.0, 2.0, -2.0, -1.0,
+ -0.5, 1.0, -1.0, -0.5
};
const int jacobianSize = spaceDim*(spaceDim-1);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2007-06-09 19:28:28 UTC (rev 7121)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2007-06-09 23:32:51 UTC (rev 7122)
@@ -40,6 +40,8 @@
{ // setUp
_data = 0;
_quadrature = 0;
+ _eqsrc = new EqKinSrc();
+ _slipfn = new BruneSlipFn();
} // setUp
// ----------------------------------------------------------------------
@@ -49,6 +51,8 @@
{ // tearDown
delete _data; _data = 0;
delete _quadrature; _quadrature = 0;
+ delete _eqsrc; _eqsrc = 0;
+ delete _slipfn; _slipfn = 0;
} // tearDown
// ----------------------------------------------------------------------
@@ -163,14 +167,18 @@
int iVertex = 0;
for (Mesh::label_sequence::iterator v_iter=vBegin;
v_iter != vEnd;
- ++v_iter) {
+ ++v_iter, ++iVertex) {
dispTpdt->updatePoint(*v_iter, &_data->fieldTpdt[iVertex*spaceDim]);
dispT->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
} // for
+ //dispT->view("DISP T");
+
// Call integrateResidual()
fault.integrateResidual(residual, &fields, mesh);
+ //residual->view("RESIDUAL");
+
// Check values
const double* valsE = _data->valsResidual;
iVertex = 0;
@@ -178,7 +186,7 @@
const double tolerance = 1.0e-06;
for (Mesh::label_sequence::iterator v_iter=vBegin;
v_iter != vEnd;
- ++v_iter) {
+ ++v_iter, ++iVertex) {
const int fiberDim = residual->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
const real_section_type::value_type* vals =
@@ -231,7 +239,7 @@
int iVertex = 0;
for (Mesh::label_sequence::iterator v_iter=vBegin;
v_iter != vEnd;
- ++v_iter) {
+ ++v_iter, ++iVertex) {
dispTpdt->updatePoint(*v_iter, &_data->fieldTpdt[iVertex*spaceDim]);
dispT->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
} // for
@@ -248,6 +256,8 @@
err = MatAssemblyEnd(jacobian, MAT_FINAL_ASSEMBLY);
CPPUNIT_ASSERT(0 == err);
+ MatView(jacobian, PETSC_VIEWER_STDOUT_WORLD);
+
const double* valsE = _data->valsJacobian;
const int nrowsE = dispT->sizeWithBC();
const int ncolsE = nrowsE;
@@ -344,7 +354,6 @@
// Check values
const double* valsE = _data->valsSlip;
- int iVertex = 0;
const int fiberDimE = spaceDim;
const double tolerance = 1.0e-06;
@@ -352,9 +361,10 @@
CPPUNIT_ASSERT(!vertices.isNull());
const Mesh::label_sequence::iterator vBegin = vertices->begin();
const Mesh::label_sequence::iterator vEnd = vertices->end();
+ int iVertex = 0;
for (Mesh::label_sequence::iterator v_iter=vBegin;
v_iter != vEnd;
- ++v_iter) {
+ ++v_iter, ++iVertex) {
const int fiberDim = disp->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
const real_section_type::value_type* vals =
@@ -379,6 +389,8 @@
CPPUNIT_ASSERT(0 != mesh);
CPPUNIT_ASSERT(0 != fault);
CPPUNIT_ASSERT(0 != _quadrature);
+ CPPUNIT_ASSERT(0 != _eqsrc);
+ CPPUNIT_ASSERT(0 != _slipfn);
try {
meshio::MeshIOAscii iohandler;
@@ -412,18 +424,16 @@
ioPeakRate.filename(_data->peakRateFilename);
dbPeakRate.ioHandler(&ioPeakRate);
- BruneSlipFn slipFn;
- slipFn.dbFinalSlip(&dbFinalSlip);
- slipFn.dbSlipTime(&dbSlipTime);
- slipFn.dbPeakRate(&dbPeakRate);
+ _slipfn->dbFinalSlip(&dbFinalSlip);
+ _slipfn->dbSlipTime(&dbSlipTime);
+ _slipfn->dbPeakRate(&dbPeakRate);
- EqKinSrc eqsrc;
- eqsrc.slipfn(&slipFn);
+ _eqsrc->slipfn(_slipfn);
fault->id(_data->id);
fault->label(_data->label);
fault->quadrature(_quadrature);
- fault->eqsrc(&eqsrc);
+ fault->eqsrc(_eqsrc);
fault->adjustTopology(*mesh);
const double upDirVals[] = { 0.0, 0.0, 1.0 };
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh 2007-06-09 19:28:28 UTC (rev 7121)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh 2007-06-09 23:32:51 UTC (rev 7122)
@@ -32,6 +32,8 @@
class FaultCohesiveKin; // USES FaultCohesiveKin
class CohesiveKinData; // HOLDSA CohesiveKinData
+ class EqKinSrc; // HOLDSA EqKinSrc
+ class BruneSlipFn; // HOLDSA BruneSlipFn
} // faults
namespace feassemble {
@@ -102,6 +104,12 @@
void _initialize(ALE::Obj<ALE::Mesh>* mesh,
FaultCohesiveKin* const fault) const;
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ EqKinSrc* _eqsrc; ///< Kinematic earthquake source
+ BruneSlipFn* _slipfn; ///< Slip time function
+
}; // class TestFaultCohesiveKin
#endif // pylith_faults_testfaultcohesivekin_hh
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc 2007-06-09 19:28:28 UTC (rev 7121)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc 2007-06-09 23:32:51 UTC (rev 7122)
@@ -27,6 +27,7 @@
void
pylith::faults::TestFaultCohesiveKinLine2::setUp(void)
{ // setUp
+ TestFaultCohesiveKin::setUp();
_data = new CohesiveKinDataLine2();
_quadrature = new feassemble::Quadrature0D();
CPPUNIT_ASSERT(0 != _quadrature);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc 2007-06-09 19:28:28 UTC (rev 7121)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc 2007-06-09 23:32:51 UTC (rev 7122)
@@ -27,6 +27,7 @@
void
pylith::faults::TestFaultCohesiveKinTri3::setUp(void)
{ // setUp
+ TestFaultCohesiveKin::setUp();
_data = new CohesiveKinDataTri3();
_quadrature = new feassemble::Quadrature1Din2D();
CPPUNIT_ASSERT(0 != _quadrature);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2007-06-09 19:28:28 UTC (rev 7121)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2007-06-09 23:32:51 UTC (rev 7122)
@@ -66,20 +66,21 @@
const char* pylith::faults::CohesiveKinDataLine2::_peakRateFilename =
"data/line2_peakrate.spatialdb";
+// Don't expect these values to be used, so just use some values.
const double pylith::faults::CohesiveKinDataLine2::_fieldTpdt[] = {
- 0.0,
- 0.0,
- 0.0,
- 0.0,
- 0.0
+ 7.1,
+ 7.2,
+ 7.3,
+ 7.4,
+ 7.5
};
const double pylith::faults::CohesiveKinDataLine2::_fieldT[] = {
- 0.0,
- 0.0,
- 0.0,
- 0.0,
- 0.0
+ 0.1,
+ 0.2,
+ 0.3,
+ 0.4,
+ 2.3 // force associated with constraint
};
const int pylith::faults::CohesiveKinDataLine2::_numConstraintVert = 1;
@@ -94,10 +95,9 @@
const double pylith::faults::CohesiveKinDataLine2::_valsResidual[] = {
0.0,
+ -2.3,
0.0,
- 0.0,
- 0.0,
- 0.0,
+ 2.3,
0.0
};
@@ -106,16 +106,15 @@
0.0,
0.0,
0.0,
- 0.0,
- 1.2
+ 1.05168389458
};
const double pylith::faults::CohesiveKinDataLine2::_valsJacobian[] = {
- 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, 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, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, +1.0,
+ 0.0, -1.0, 0.0, +1.0, 0.0,
};
pylith::faults::CohesiveKinDataLine2::CohesiveKinDataLine2(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2007-06-09 19:28:28 UTC (rev 7121)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2007-06-09 23:32:51 UTC (rev 7122)
@@ -93,32 +93,32 @@
"data/tri3_peakrate.spatialdb";
const double pylith::faults::CohesiveKinDataTri3::_fieldTpdt[] = {
- 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,
+ 8.1, 9.1,
+ 8.2, 9.2,
+ 8.3, 9.3,
+ 8.4, 9.4,
+ 8.5, 9.5,
+ 8.6, 9.6,
+ 8.7, 9.7,
+ 8.8, 9.8,
};
const double pylith::faults::CohesiveKinDataTri3::_fieldT[] = {
- 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,
+ 6.1, 7.1, // 2
+ 6.2, 7.2, // 3
+ 6.3, 7.3, // 4
+ 6.4, 7.4, // 5
+ 6.5, 7.5, // 6
+ 1.4, 1.5, // 7 (constraint force)
+ 6.7, 7.7, // 8
+ 2.6, 2.7, // 9 (constraint force)
};
const int pylith::faults::CohesiveKinDataTri3::_numConstraintVert = 2;
const double pylith::faults::CohesiveKinDataTri3::_orientation[] = {
- 0.0, -1.0, -1.0, 0.0,
- 0.0, -1.0, -1.0, 0.0
+ 0.0, -1.0, +1.0, 0.0,
+ 0.0, -1.0, +1.0, 0.0
};
const int pylith::faults::CohesiveKinDataTri3::_constraintVertices[] = {
@@ -126,25 +126,25 @@
};
const double pylith::faults::CohesiveKinDataTri3::_valsResidual[] = {
- 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, // 2
+ -1.4, -1.5, // 3
+ -2.6, -2.7, // 4
+ 0.0, 0.0, // 5
+ +1.4, +1.5, // 6
+ 0.0, 0.0, // 7
+ +2.6, +2.7, // 8
+ 0.0, 0.0 // 9
};
const double pylith::faults::CohesiveKinDataTri3::_valsSlip[] = {
- 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,
+ 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,
};
const double pylith::faults::CohesiveKinDataTri3::_valsJacobian[] = {
More information about the cig-commits
mailing list