[cig-commits] r16443 - in short/3D/PyLith/trunk: libsrc/faults libsrc/feassemble libsrc/problems libsrc/topology playpen/friction/bar_shearwave/quad4 playpen/friction/twohex8 pylith/problems unittests/libtests/bc unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Tue Mar 23 15:52:58 PDT 2010
Author: brad
Date: 2010-03-23 15:52:56 -0700 (Tue, 23 Mar 2010)
New Revision: 16443
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh
short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/pylithapp.cfg
short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/shear-sliding.cfg
short/3D/PyLith/trunk/playpen/friction/twohex8/opening.cfg
short/3D/PyLith/trunk/pylith/problems/Explicit.py
short/3D/PyLith/trunk/pylith/problems/Formulation.py
short/3D/PyLith/trunk/pylith/problems/Implicit.py
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
Log:
Implemented sensitivity solve (perform linear solve to determine change in slip associated with change in Lagrange multipliers) using vertices associated with fault. Refactor friction methods specific to spatial dimension. Still need to update unit tests.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2010-03-23 22:52:56 UTC (rev 16443)
@@ -116,7 +116,11 @@
bool _useLagrangeConstraints; ///< True if uses Lagrange multipliers.
- // PRIVATE MEMBERS ////////////////////////////////////////////////////
+ /// Map label of cohesive cell to label of fault cell.
+ std::map<topology::Mesh::SieveMesh::point_type,
+ topology::SubMesh::SieveMesh::point_type> _cohesiveToFault;
+
+// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
/// If true, use fault mesh to define fault; otherwise, use group of
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-03-23 22:52:56 UTC (rev 16443)
@@ -26,6 +26,7 @@
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/friction/FrictionModel.hh" // USES FrictionModel
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
+#include "pylith/problems/SolverLinear.hh" // USES SolverLinear
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
@@ -53,9 +54,10 @@
// Default constructor.
pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(void) :
_dbInitialTract(0),
- _friction(0)
+ _friction(0),
+ _jacobian(0),
+ _ksp(0)
{ // constructor
- _needJacobianDiag = true;
_needVelocity = true;
} // constructor
@@ -74,6 +76,12 @@
_dbInitialTract = 0; // :TODO: Use shared pointer
_friction = 0; // :TODO: Use shared pointer
+
+ delete _jacobian; _jacobian = 0;
+ if (0 != _ksp) {
+ PetscErrorCode err = KSPDestroy(_ksp); _ksp = 0;
+ CHECK_PETSC_ERROR(err);
+ } // if
} // deallocate
// ----------------------------------------------------------------------
@@ -316,11 +324,12 @@
const double t,
const topology::Jacobian& jacobian)
{ // constrainSolnSpace
- /// Member prototype for _constrainSolnSpaceLumpedXD()
+ /// Member prototype for _constrainSolnSpaceXD()
typedef void (pylith::faults::FaultCohesiveDyn::*constrainSolnSpace_fn_type)
- (double_array*, double_array*, const double_array&,
- const double_array&, const double_array&,
- const double_array&, const double_array&,
+ (double_array*,
+ const double_array&,
+ const double_array&,
+ const double_array&,
const double);
assert(0 != fields);
@@ -329,6 +338,7 @@
assert(0 != _friction);
_updateSlipRate(*fields);
+ _sensitivitySetup(jacobian);
const int spaceDim = _quadrature->spaceDim();
@@ -337,10 +347,10 @@
double_array tractionTpdtVertex(spaceDim);
double_array slipTpdtVertex(spaceDim);
double_array lagrangeTpdtVertex(spaceDim);
- double_array dLagrangeTpdtVertex(spaceDim);
// Get sections
double_array slipVertex(spaceDim);
+ double_array dSlipVertex(spaceDim);
const ALE::Obj<RealSection>& slipSection = _fields->get("slip").section();
assert(!slipSection.isNull());
@@ -358,33 +368,36 @@
assert(!orientationSection.isNull());
double_array lagrangeTVertex(spaceDim);
+ double_array dispTVertexN(spaceDim);
+ double_array dispTVertexP(spaceDim);
const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
double_array lagrangeTIncrVertex(spaceDim);
+ double_array dispTIncrVertexN(spaceDim);
+ double_array dispTIncrVertexP(spaceDim);
const ALE::Obj<RealSection>& dispTIncrSection =
fields->get("dispIncr(t->t+dt)").section();
assert(!dispTIncrSection.isNull());
- double_array jacobianVertexN(spaceDim);
- double_array jacobianVertexP(spaceDim);
- const ALE::Obj<RealSection>& jacobianDiagSection =
- fields->get("Jacobian diagonal").section();
- assert(!jacobianDiagSection.isNull());
+ double_array dLagrangeTpdtVertex(spaceDim);
+ const ALE::Obj<RealSection>& dLagrangeTpdtSection =
+ _fields->get("sensitivity dLagrange").section();
+ assert(!dLagrangeTpdtSection.isNull());
constrainSolnSpace_fn_type constrainSolnSpaceFn;
switch (spaceDim) { // switch
case 1:
constrainSolnSpaceFn =
- &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped1D;
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D;
break;
case 2:
constrainSolnSpaceFn =
- &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped2D;
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D;
break;
case 3:
constrainSolnSpaceFn =
- &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped3D;
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D;
break;
default :
assert(0);
@@ -404,8 +417,6 @@
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
- const int v_negative = _cohesiveVertices[iVertex].negative;
- const int v_positive = _cohesiveVertices[iVertex].positive;
// Get slip
slipSection->restrictPoint(v_fault, &slipVertex[0], slipVertex.size());
@@ -419,17 +430,6 @@
assert(0 != areaVertex);
assert(1 == areaSection->getFiberDimension(v_fault));
- // Get fault orientation
- orientationSection->restrictPoint(v_fault, &orientationVertex[0],
- orientationVertex.size());
-
- // Get diagonal of Jacobian at conventional vertices i and j
- // associated with Lagrange vertex k
- jacobianDiagSection->restrictPoint(v_negative, &jacobianVertexN[0],
- jacobianVertexN.size());
- jacobianDiagSection->restrictPoint(v_positive, &jacobianVertexP[0],
- jacobianVertexP.size());
-
// Get Lagrange multiplier values from disp(t), and dispIncr(t->t+dt)
dispTSection->restrictPoint(v_lagrange, &lagrangeTVertex[0],
lagrangeTVertex.size());
@@ -454,25 +454,85 @@
// friction.
CALL_MEMBER_FN(*this,
constrainSolnSpaceFn)(&dLagrangeTpdtVertex,
- &slipVertex, slipRateVertex,
- tractionTpdtVertex, orientationVertex,
- jacobianVertexN, jacobianVertexP,
- *areaVertex);
+ slipVertex, slipRateVertex,
+ tractionTpdtVertex, *areaVertex);
- // Update Lagrange multiplier values.
- lagrangeTIncrVertex += dLagrangeTpdtVertex;
- assert(lagrangeTIncrVertex.size() ==
- dispTIncrSection->getFiberDimension(v_lagrange));
- dispTIncrSection->updatePoint(v_lagrange, &lagrangeTIncrVertex[0]);
+ assert(dLagrangeTpdtVertex.size() ==
+ dLagrangeTpdtSection->getFiberDimension(v_fault));
+ dLagrangeTpdtSection->updatePoint(v_fault, &dLagrangeTpdtVertex[0]);
+ } // for
- // Update the slip estimate based on adjustment to the Lagrange
- // multiplier values.
- assert(slipVertex.size() ==
+ // Solve sensitivity problem for negative side of the fault.
+ bool negativeSide = true;
+ _sensitivityUpdateJacobian(negativeSide, jacobian, *fields);
+ _sensitivityReformResidual(negativeSide);
+ _sensitivitySolve();
+ _sensitivityUpdateSoln(negativeSide);
+
+ // Solve sensitivity problem for positive side of the fault.
+ negativeSide = false;
+ _sensitivityUpdateJacobian(negativeSide, jacobian, *fields);
+ _sensitivityReformResidual(negativeSide);
+ _sensitivitySolve();
+ _sensitivityUpdateSoln(negativeSide);
+
+ // Update slip field based on solution of sensitivity problem and
+ // increment in Lagrange multipliers.
+ const ALE::Obj<RealSection>& dispRelSection =
+ _fields->get("sensitivity dispRel").section();
+ for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+ const int v_fault = _cohesiveVertices[iVertex].fault;
+ const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
+ const int v_negative = _cohesiveVertices[iVertex].negative;
+ const int v_positive = _cohesiveVertices[iVertex].positive;
+
+ // Get fault orientation
+ orientationSection->restrictPoint(v_fault, &orientationVertex[0],
+ orientationVertex.size());
+
+ // Get change in relative displacement.
+ const double* dispRelVertex = dispRelSection->restrictPoint(v_fault);
+ assert(0 != dispRelVertex);
+ assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
+
+ // Get Lagrange multiplier at time t
+ dispTSection->restrictPoint(v_lagrange, &lagrangeTVertex[0],
+ lagrangeTVertex.size());
+
+ // Get Lagrange multiplier increment at time t
+ dispTIncrSection->restrictPoint(v_lagrange, &lagrangeTIncrVertex[0],
+ lagrangeTIncrVertex.size());
+
+ // Get change in Lagrange multiplier.
+ dLagrangeTpdtSection->restrictPoint(v_fault, &dLagrangeTpdtVertex[0],
+ dLagrangeTpdtVertex.size());
+ // Compute change in slip.
+ dSlipVertex = 0.0;
+ for (int iDim = 0; iDim < spaceDim; ++iDim)
+ for (int kDim = 0; kDim < spaceDim; ++kDim)
+ dSlipVertex[iDim] +=
+ orientationVertex[iDim*spaceDim+kDim] * dispRelVertex[kDim];
+
+ // Set fault opening to zero if fault is under compression.
+ const int indexN = spaceDim - 1;
+ const double lagrangeTpdtNormal = lagrangeTVertex[indexN] +
+ lagrangeTIncrVertex[indexN] + dLagrangeTpdtVertex[indexN];
+ if (lagrangeTpdtNormal < 0.0)
+ dSlipVertex[indexN] = 0.0;
+
+ // Set change in slip.
+ assert(dSlipVertex.size() ==
slipSection->getFiberDimension(v_fault));
- slipSection->updatePoint(v_fault, &slipVertex[0]);
- } // if
+ slipSection->updateAddPoint(v_fault, &dSlipVertex[0]);
+
+ // Update Lagrange multiplier increment.
+ assert(dLagrangeTpdtVertex.size() ==
+ dispTIncrSection->getFiberDimension(v_lagrange));
+ dispTIncrSection->updateAddPoint(v_lagrange, &dLagrangeTpdtVertex[0]);
+ } // for
#if 0 // DEBUGGING
+ dLagrangeTpdtSection->view("AFTER dLagrange");
dispTIncrSection->view("AFTER DISP INCR (t->t+dt)");
slipSection->view("AFTER SLIP");
#endif
@@ -486,19 +546,28 @@
topology::SolutionFields* const fields,
const topology::Field<topology::Mesh>& jacobian)
{ // adjustSolnLumped
- /// Member prototype for _constrainSolnSpaceLumpedXD()
+ /// Member prototype for _constrainSolnSpaceXD()
typedef void (pylith::faults::FaultCohesiveDyn::*constrainSolnSpace_fn_type)
- (double_array*, double_array*, const double_array&,
- const double_array&, const double_array&,
- const double_array&, const double_array&,
+ (double_array*,
+ const double_array&,
+ const double_array&,
+ const double_array&,
const double);
+ /// Member prototype for _sensitivitySolveLumpedXD()
+ typedef void (pylith::faults::FaultCohesiveDyn::*sensitivitySolveLumped_fn_type)
+ (double_array*,
+ const double_array&,
+ const double_array&,
+ const double_array&,
+ const double_array&);
+
/// Member prototype for _adjustSolnLumpedXD()
typedef void (pylith::faults::FaultCohesiveDyn::*adjustSolnLumped_fn_type)
(double_array*, double_array*, double_array*,
- const double_array&, const double_array&,
- const double_array&, const double_array&,
- const double_array&, const double_array&,
+ const double_array&, const double_array&,
+ const double_array&, const double_array&,
+ const double_array&, const double_array&,
const double_array&, const double_array&);
@@ -581,24 +650,31 @@
adjustSolnLumped_fn_type adjustSolnLumpedFn;
constrainSolnSpace_fn_type constrainSolnSpaceFn;
+ sensitivitySolveLumped_fn_type sensitivitySolveLumpedFn;
switch (spaceDim) { // switch
case 1:
adjustSolnLumpedFn =
&pylith::faults::FaultCohesiveDyn::_adjustSolnLumped1D;
constrainSolnSpaceFn =
- &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped1D;
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D;
+ sensitivitySolveLumpedFn =
+ &pylith::faults::FaultCohesiveDyn::_sensitivitySolveLumped1D;
break;
case 2:
adjustSolnLumpedFn =
&pylith::faults::FaultCohesiveDyn::_adjustSolnLumped2D;
constrainSolnSpaceFn =
- &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped2D;
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D;
+ sensitivitySolveLumpedFn =
+ &pylith::faults::FaultCohesiveDyn::_sensitivitySolveLumped2D;
break;
case 3:
adjustSolnLumpedFn =
&pylith::faults::FaultCohesiveDyn::_adjustSolnLumped3D;
constrainSolnSpaceFn =
- &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped3D;
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D;
+ sensitivitySolveLumpedFn =
+ &pylith::faults::FaultCohesiveDyn::_sensitivitySolveLumped3D;
break;
default :
assert(0);
@@ -699,13 +775,16 @@
CALL_MEMBER_FN(*this,
constrainSolnSpaceFn)(&dLagrangeTpdtVertex,
- &slipVertex, slipRateVertex,
- tractionTpdtVertex, orientationVertex,
- jacobianVertexN, jacobianVertexP,
- *areaVertex);
+ slipVertex, slipRateVertex,
+ tractionTpdtVertex, *areaVertex);
+ CALL_MEMBER_FN(*this,
+ sensitivitySolveLumpedFn)(&slipVertex,
+ dLagrangeTpdtVertex, orientationVertex,
+ jacobianVertexN, jacobianVertexP);
lagrangeTIncrVertex += dLagrangeTpdtVertex;
+ // :TODO: Refactor this into sensitivitySolveLumpedXD().
switch (spaceDim) { // switch
case 1: {
assert(jacobianVertexN[0] > 0.0);
@@ -1340,13 +1419,13 @@
for (int iDim = 0; iDim < spaceDim; ++iDim)
for (int kDim = 0; kDim < spaceDim; ++kDim)
slipRateVertex[iDim] +=
- velocityVertexN[kDim] * -orientationVertex[kDim*spaceDim+iDim];
+ velocityVertexN[kDim] * -orientationVertex[iDim*spaceDim+kDim];
// Velocity for positive vertex.
for (int iDim = 0; iDim < spaceDim; ++iDim)
for (int kDim = 0; kDim < spaceDim; ++kDim)
slipRateVertex[iDim] +=
- velocityVertexP[kDim] * +orientationVertex[kDim*spaceDim+iDim];
+ velocityVertexP[kDim] * +orientationVertex[iDim*spaceDim+kDim];
// Update slip rate field.
assert(slipRateVertex.size() == slipRateSection->getFiberDimension(v_fault));
@@ -1357,94 +1436,323 @@
} // _updateSlipRate
// ----------------------------------------------------------------------
-// Update Jacobian blocks associated with DOF of vertices on negative
-// and positive sides of the fault associated with Lagrange vertex k.
+// Setup sensitivity problem to compute change in slip given change in Lagrange multipliers.
void
-pylith::faults::FaultCohesiveDyn::_updateJacobianBlocks(
- const topology::Jacobian& jacobian,
- const topology::SolutionFields& fields)
-{ // _updateJacobianBlocks
+pylith::faults::FaultCohesiveDyn::_sensitivitySetup(const topology::Jacobian& jacobian)
+{ // _sensitivitySetup
assert(0 != _fields);
assert(0 != _quadrature);
const int spaceDim = _quadrature->spaceDim();
- if (!_fields->hasField("Jacobian blocks")) {
- // Create field for entries of Jacobian at conventional vertices on
- // negative and positive sides of the fault associated with Lagrange
- // vertex k.
- _fields->add("Jacobian blocks", "jacobian_blocks");
- topology::Field<topology::SubMesh>& jacobianBlocks =
- _fields->get("Jacobian diagonal");
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- jacobianBlocks.newSection(slip, 2*spaceDim*spaceDim);
- jacobianBlocks.allocate();
- jacobianBlocks.vectorFieldType(topology::FieldBase::OTHER);
+ // Setup fields involved in sensitivity solve.
+ if (!_fields->hasField("sensitivity solution")) {
+ _fields->add("sensitivity solution", "sensitivity_soln");
+ topology::Field<topology::SubMesh>& solution =
+ _fields->get("sensitivity solution");
+ const topology::Field<topology::SubMesh>& slip =
+ _fields->get("slip");
+ solution.cloneSection(slip);
+ solution.createVector();
+ solution.createScatter();
} // if
-
- const ALE::Obj<RealSection>& jacobianBlocksSection =
- _fields->get("Jacobian diagonal").section();
+ const topology::Field<topology::SubMesh>& solution =
+ _fields->get("sensitivity solution");
+ if (!_fields->hasField("sensitivity residual")) {
+ _fields->add("sensitivity residual", "sensitivity_residual");
+ topology::Field<topology::SubMesh>& residual =
+ _fields->get("sensitivity residual");
+ residual.cloneSection(solution);
+ residual.createVector();
+ residual.createScatter();
+ } // if
+
+ if (!_fields->hasField("sensitivity dispRel")) {
+ _fields->add("sensitivity dispRel", "sensitivity_disprel");
+ topology::Field<topology::SubMesh>& dispRel =
+ _fields->get("sensitivity dispRel");
+ dispRel.cloneSection(solution);
+ } // if
+ topology::Field<topology::SubMesh>& dispRel =
+ _fields->get("sensitivity dispRel");
+ dispRel.zero();
+
+ if (!_fields->hasField("sensitivity dLagrange")) {
+ _fields->add("sensitivity dLagrange", "sensitivity_dlagrange");
+ topology::Field<topology::SubMesh>& dLagrange =
+ _fields->get("sensitivity dLagrange");
+ dLagrange.cloneSection(solution);
+ } // if
+ topology::Field<topology::SubMesh>& dLagrange =
+ _fields->get("sensitivity dLagrange");
+ dLagrange.zero();
+
+ // Setup Jacobian sparse matrix for sensitivity solve.
+ if (0 == _jacobian)
+ _jacobian = new topology::Jacobian(solution, jacobian.matrixType());
+ assert(0 != _jacobian);
+ _jacobian->zero();
+
+ // Setup PETSc KSP linear solver.
+ if (0 == _ksp) {
+ PetscErrorCode err = 0;
+ if (0 != _ksp) {
+ err = KSPDestroy(_ksp); _ksp = 0;
+ CHECK_PETSC_ERROR(err);
+ } // if
+ err = KSPCreate(_faultMesh->comm(), &_ksp); CHECK_PETSC_ERROR(err);
+ err = KSPSetInitialGuessNonzero(_ksp, PETSC_FALSE); CHECK_PETSC_ERROR(err);
+
+ PC pc;
+ err = KSPGetPC(_ksp, &pc); CHECK_PETSC_ERROR(err);
+ err = PCSetType(pc, PCJACOBI); CHECK_PETSC_ERROR(err);
+ err = KSPSetType(_ksp, KSPGMRES);
+
+ err = KSPSetFromOptions(_ksp); CHECK_PETSC_ERROR(err);
+ err = KSPAppendOptionsPrefix(_ksp, "friction");
+ } // if
+} // _sensitivitySetup
+
+// ----------------------------------------------------------------------
+// Update the Jacobian values for the sensitivity solve.
+void
+pylith::faults::FaultCohesiveDyn::_sensitivityUpdateJacobian(const bool negativeSide,
+ const topology::Jacobian& jacobian,
+ const topology::SolutionFields& fields)
+{ // _sensitivityUpdateJacobian
+ assert(0 != _quadrature);
+ assert(0 != _fields);
+
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int subnrows = numBasis*spaceDim;
+ const int submatrixSize = subnrows * subnrows;
+ const int nrows = 3*subnrows;
+ const int matrixSize = nrows * nrows;
+
+ // Get solution field
+ const topology::Field<topology::Mesh>& solutionDomain = fields.solution();
+ const ALE::Obj<RealSection>& solutionDomainSection = solutionDomain.section();
+ assert(!solutionDomainSection.isNull());
+
+ // Get cohesive cells
const ALE::Obj<SieveMesh>& sieveMesh = fields.mesh().sieveMesh();
assert(!sieveMesh.isNull());
- const ALE::Obj<RealSection>& solutionSection = fields.solution().section();
- assert(!solutionSection.isNull());
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
- solutionSection);
- assert(!globalOrder.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
+ sieveMesh->getLabelStratum("material-id", id());
+ assert(!cellsCohesive.isNull());
+ const SieveMesh::label_sequence::iterator cellsCohesiveBegin =
+ cellsCohesive->begin();
+ const SieveMesh::label_sequence::iterator cellsCohesiveEnd =
+ cellsCohesive->end();
- const PetscMat jacobianMatrix = jacobian.matrix();
- assert(0 != jacobianMatrix);
+ // Visitor for Jacobian matrix associated with domain.
+ double_array jacobianDomainCell(matrixSize);
+ const PetscMat jacobianDomainMatrix = jacobian.matrix();
+ assert(0 != jacobianDomainMatrix);
+ const ALE::Obj<SieveMesh::order_type>& globalOrderDomain =
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solutionDomainSection);
+ assert(!globalOrderDomain.isNull());
+ // We would need to request unique points here if we had an interpolated mesh
+ topology::Mesh::IndicesVisitor jacobianDomainVisitor(*solutionDomainSection,
+ *globalOrderDomain,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
- // Arrays for values (extracting block diagonal so rows and cols are
- // the same).
- const int nvalues = 2*spaceDim*spaceDim;
- const int nrows = 2*spaceDim;
- double_array jacobianValues(nvalues);
- int_array rows(nrows);
+ // Get fault Sieve mesh
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ // Get sensitivity solution field
+ const ALE::Obj<RealSection>& solutionFaultSection =
+ _fields->get("sensitivity solution").section();
+ assert(!solutionFaultSection.isNull());
+
+ // Visitor for Jacobian matrix associated with fault.
+ double_array jacobianFaultCell(submatrixSize);
+ assert(0 != _jacobian);
+ const PetscMat jacobianFaultMatrix = _jacobian->matrix();
+ assert(0 != jacobianFaultMatrix);
+ const ALE::Obj<SieveSubMesh::order_type>& globalOrderFault =
+ faultSieveMesh->getFactory()->getGlobalOrder(faultSieveMesh, "default", solutionFaultSection);
+ assert(!globalOrderFault.isNull());
+ // We would need to request unique points here if we had an interpolated mesh
+ topology::SubMesh::IndicesVisitor jacobianFaultVisitor(*solutionFaultSection,
+ *globalOrderFault,
+ (int) pow(faultSieveMesh->getSieve()->getMaxConeSize(),
+ faultSieveMesh->depth())*spaceDim);
+
+ const int indexOffset = (negativeSide) ? 0 : 3*submatrixSize + subnrows;
+
+ for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
+ c_iter != cellsCohesiveEnd;
+ ++c_iter) {
+ const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
+ jacobianDomainCell = 0.0;
+ jacobianFaultCell = 0.0;
+
+ // Get cell contribution in PETSc Matrix
+ jacobianDomainVisitor.clear();
+#if 0
+
+ PetscErrorCode err = restrictOperator(jacobianDomainMatrix, *sieveMesh->getSieve(),
+ jacobianDomainVisitor, *c_iter,
+ &jacobianDomainCell[0]);
+ CHECK_PETSC_ERROR_MSG(err, "Restrict from PETSc Mat failed.");
+#else
+ ALE::ISieveTraversal<SieveMesh::sieve_type>::orientedClosure(*sieveMesh->getSieve(),
+ *c_iter,
+ jacobianDomainVisitor);
+ const int* indices = jacobianDomainVisitor.getValues();
+ const int numIndices = jacobianDomainVisitor.getSize();
+ PetscErrorCode err = MatGetValues(jacobianDomainMatrix,
+ numIndices, indices,
+ numIndices, indices,
+ &jacobianDomainCell[0]);
+ CHECK_PETSC_ERROR_MSG(err, "Restrict from PETSc Mat failed.");
+#endif
+
+ for (int iRow=0, i=0; iRow < subnrows; ++iRow) {
+ const int indexR = indexOffset + iRow*nrows;
+ for (int iCol=0; iCol < subnrows; ++iCol, ++i)
+ jacobianFaultCell[i] = jacobianDomainCell[indexR+iCol];
+ } // for
+
+ // Insert cell contribution into PETSc Matrix
+ jacobianFaultVisitor.clear();
+ err = updateOperator(jacobianFaultMatrix, *faultSieveMesh->getSieve(),
+ jacobianFaultVisitor, c_fault,
+ &jacobianFaultCell[0], INSERT_VALUES);
+ CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
+ } // for
+
+ _jacobian->assemble("final_assembly");
+} // _sensitivityUpdateJacobian
+
+// ----------------------------------------------------------------------
+// Reform residual for sensitivity problem.
+void
+pylith::faults::FaultCohesiveDyn::_sensitivityReformResidual(const bool negativeSide)
+{ // _sensitivityReformResidual
+ assert(0 != _fields);
+ assert(0 != _quadrature);
+
+ const int spaceDim = _quadrature->spaceDim();
+
+ // Compute residual -C^T dLagrange
+ double_array residualVertex(spaceDim);
+ topology::Field<topology::SubMesh>& residual =
+ _fields->get("sensitivity residual");
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ residual.zero();
+
+ const ALE::Obj<RealSection>& dLagrangeSection =
+ _fields->get("sensitivity dLagrange").section();
+
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
+ assert(!orientationSection.isNull());
+
+ const double sign = (negativeSide) ? -1.0 : 1.0;
+
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_fault = _cohesiveVertices[iVertex].fault;
- const int v_negative = _cohesiveVertices[iVertex].negative;
- const int v_positive = _cohesiveVertices[iVertex].positive;
- // Set rows/cols using globalOrder
- const int indexN = globalOrder->getIndex(v_negative);
- const int indexP = globalOrder->getIndex(v_positive);
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- rows[0*spaceDim+iDim] = indexN + iDim;
- rows[1*spaceDim+iDim] = indexP + iDim;
- } // for
+ const double* dLagrangeVertex = dLagrangeSection->restrictPoint(v_fault);
+ assert(0 != dLagrangeVertex);
+ assert(spaceDim == dLagrangeSection->getFiberDimension(v_fault));
- // Get values
- PetscErrorCode err = MatGetValues(jacobianMatrix,
- nrows, &rows[0],
- nrows, &rows[0],
- &jacobianValues[0]);
-
- // Store blocks
- assert(jacobianValues.size() ==
- jacobianBlocksSection->getFiberDimension(v_fault));
- jacobianBlocksSection->updatePoint(v_fault, &jacobianValues[0]);
+ const double* orientationVertex = orientationSection->restrictPoint(v_fault);
+ assert(0 != orientationVertex);
+ assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+
+ residualVertex = 0.0;
+ for (int iDim = 0; iDim < spaceDim; ++iDim)
+ for (int kDim = 0; kDim < spaceDim; ++kDim)
+ residualVertex[iDim] +=
+ sign * dLagrangeVertex[kDim] * -orientationVertex[kDim*spaceDim+iDim];
+
+ assert(residualVertex.size() == residualSection->getFiberDimension(v_fault));
+ residualSection->updatePoint(v_fault, &residualVertex[0]);
} // for
-} // _updateJacobianBlocks
+ PetscLogFlops(numVertices*spaceDim*spaceDim*4);
+} // _sensitivityReformResidual
+
// ----------------------------------------------------------------------
-// Constrain solution space in 1-D.
+// Solve sensitivity problem.
void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped1D(
- double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& sliprate,
- const double_array& tractionTpdt,
+pylith::faults::FaultCohesiveDyn::_sensitivitySolve(void)
+{ // _sensitivitySolve
+ assert(0 != _fields);
+ assert(0 != _jacobian);
+ assert(0 != _ksp);
+
+ const topology::Field<topology::SubMesh>& residual =
+ _fields->get("sensitivity residual");
+ const topology::Field<topology::SubMesh>& solution =
+ _fields->get("sensitivity solution");
+
+ // Update PetscVector view of field.
+ residual.scatterSectionToVector();
+
+ PetscErrorCode err = 0;
+ const PetscMat jacobianMat = _jacobian->matrix();
+ err = KSPSetOperators(_ksp, jacobianMat, jacobianMat,
+ DIFFERENT_NONZERO_PATTERN); CHECK_PETSC_ERROR(err);
+
+ const PetscVec residualVec = residual.vector();
+ const PetscVec solutionVec = solution.vector();
+ err = KSPSolve(_ksp, residualVec, solutionVec); CHECK_PETSC_ERROR(err);
+
+ // Update section view of field.
+ solution.scatterVectorToSection();
+} // _sensitivitySolve
+
+// ----------------------------------------------------------------------
+// Update the relative displacement field values based on the
+// sensitivity solve.
+void
+pylith::faults::FaultCohesiveDyn::_sensitivityUpdateSoln(const bool negativeSide)
+{ // _sensitivityUpdateSoln
+ assert(0 != _fields);
+ assert(0 != _quadrature);
+
+ const int spaceDim = _quadrature->spaceDim();
+
+ double_array dispVertex(spaceDim);
+ const ALE::Obj<RealSection>& solutionSection =
+ _fields->get("sensitivity solution").section();
+ const ALE::Obj<RealSection>& dispRelSection =
+ _fields->get("sensitivity dispRel").section();
+
+ const double sign = (negativeSide) ? -1.0 : 1.0;
+
+ const int numVertices = _cohesiveVertices.size();
+ for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+ const int v_fault = _cohesiveVertices[iVertex].fault;
+
+ solutionSection->restrictPoint(v_fault, &dispVertex[0], dispVertex.size());
+
+ dispVertex *= sign;
+
+ assert(dispVertex.size() == dispRelSection->getFiberDimension(v_fault));
+ dispRelSection->updateAddPoint(v_fault, &dispVertex[0]);
+ } // for
+} // _sensitivityUpdateSoln
+
+// ----------------------------------------------------------------------
+// Solve slip/Lagrange multiplier sensitivity problem for case of lumped Jacobian in 1-D.
+void
+pylith::faults::FaultCohesiveDyn::_sensitivitySolveLumped1D(
+ double_array* slip,
+ const double_array& dLagrangeTpdt,
const double_array& orientation,
const double_array& jacobianN,
- const double_array& jacobianP,
- const double area)
-{ // constrainSolnSpace1D
- assert(0 != dLagrangeTpdt);
+ const double_array& jacobianP)
+{ // _sensitivitySolveLumped1D
assert(0 != slip);
// Sensitivity of slip to changes in the Lagrange multipliers
@@ -1452,40 +1760,27 @@
const double Aixjx = 1.0 / jacobianN[0] + 1.0
/ jacobianP[0];
const double Spp = 1.0;
-
- if (tractionTpdt[0] < 0) {
- // if compression, then no changes to solution
- } else {
- // if tension, then traction is zero.
-
- // Update slip based on value required to stick versus
- // zero traction
- const double dlp = tractionTpdt[0] * area;
- (*dLagrangeTpdt)[0] = -dlp;
- (*slip)[0] += Spp * dlp;
- } // else
- PetscLogFlops(6);
-} // constrainSolnSpace1D
+ const double dlp = dLagrangeTpdt[0];
+ (*slip)[0] -= Spp * dlp;
+ std::cout << "Slip: (" << (*slip)[0] << ")" << std::endl;
+
+ PetscLogFlops(2);
+} // _sensitivitySolveLumped1D
+
// ----------------------------------------------------------------------
-// Constrain solution space in 2-D.
+// Solve slip/Lagrange multiplier sensitivity problem for case of lumped Jacobian in 2-D.
void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped2D(
- double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& slipRate,
- const double_array& tractionTpdt,
+pylith::faults::FaultCohesiveDyn::_sensitivitySolveLumped2D(
+ double_array* slip,
+ const double_array& dLagrangeTpdt,
const double_array& orientation,
const double_array& jacobianN,
- const double_array& jacobianP,
- const double area)
-{ // constrainSolnSpace2D
- assert(0 != dLagrangeTpdt);
+ const double_array& jacobianP)
+{ // _sensitivitySolveLumped2D
assert(0 != slip);
- std::cout << "Normal traction:" << tractionTpdt[1] << std::endl;
-
// Sensitivity of slip to changes in the Lagrange multipliers
// Aixjx = 1.0/Aix + 1.0/Ajx
assert(jacobianN[0] > 0.0);
@@ -1506,71 +1801,29 @@
const double Spp = Cpx * Cpx * Aixjx + Cpy * Cpy * Aiyjy;
const double Spq = Cpx * Cqx * Aixjx + Cpy * Cqy * Aiyjy;
const double Sqq = Cqx * Cqx * Aixjx + Cqy * Cqy * Aiyjy;
-
- const double slipMag = fabs((*slip)[0]);
- const double slipRateMag = fabs(slipRate[0]);
- const double tractionNormal = tractionTpdt[1];
- const double tractionShearMag = fabs(tractionTpdt[0]);
+ const double dlp = dLagrangeTpdt[0];
+ const double dlq = dLagrangeTpdt[1];
+ (*slip)[0] -= Spp * dlp + Spq * dlq;
+ (*slip)[1] -= Spq * dlp + Sqq * dlq;
- if (tractionNormal < 0 && 0.0 == (*slip)[1]) {
- // if in compression and no opening
- std::cout << "FAULT IN COMPRESSION" << std::endl;
- const double frictionStress = _friction->calcFriction(slipMag, slipRateMag,
- tractionNormal);
- std::cout << "frictionStress: " << frictionStress << std::endl;
- if (tractionShearMag > frictionStress ||
- (tractionShearMag < frictionStress && slipMag > 0.0)) {
- // traction is limited by friction, so have sliding
- std::cout << "LIMIT TRACTION, HAVE SLIDING" << std::endl;
-
- // Update slip based on value required to stick versus friction
- const double dlp = (tractionShearMag - frictionStress) * area *
- tractionTpdt[0] / tractionShearMag;
- (*dLagrangeTpdt)[0] = -dlp;
- (*slip)[0] += Spp * dlp;
- std::cout << "Estimated slip: " << (*slip)[0] << std::endl;
- } else {
- // else friction exceeds value necessary, so stick
- std::cout << "STICK" << std::endl;
- // no changes to solution
- } // if/else
- } else {
- // if in tension, then traction is zero.
- std::cout << "FAULT IN TENSION" << std::endl;
-
- // Update slip based on value required to stick versus
- // zero traction
- const double dlp = tractionTpdt[0] * area;
- const double dlq = tractionTpdt[1] * area;
+ std::cout << "Slip: (" << (*slip)[0] << ", " << (*slip)[1] << ")" << std::endl;
- (*dLagrangeTpdt)[0] = -dlp;
- (*dLagrangeTpdt)[1] = -dlq;
- (*slip)[0] += Spp * dlp + Spq * dlq;
- (*slip)[1] += Spq * dlp + Sqq * dlq;
- } // else
+ PetscLogFlops(8);
+} // _sensitivitySolveLumped2D
- PetscLogFlops(0); // :TODO: Fix this
-} // constrainSolnSpace2D
-
// ----------------------------------------------------------------------
-// Constrain solution space in 3-D.
+// Solve slip/Lagrange multiplier sensitivity problem for case of lumped Jacobian in 3-D.
void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped3D(
- double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& slipRate,
- const double_array& tractionTpdt,
+pylith::faults::FaultCohesiveDyn::_sensitivitySolveLumped3D(
+ double_array* slip,
+ const double_array& dLagrangeTpdt,
const double_array& orientation,
const double_array& jacobianN,
- const double_array& jacobianP,
- const double area)
-{ // constrainSolnSpace3D
- assert(0 != dLagrangeTpdt);
+ const double_array& jacobianP)
+{ // _sensitivitySolveLumped3D
assert(0 != slip);
- std::cout << "Normal traction:" << tractionTpdt[2] << std::endl;
-
// Sensitivity of slip to changes in the Lagrange multipliers
// Aixjx = 1.0/Aix + 1.0/Ajx
assert(jacobianN[0] > 0.0);
@@ -1595,47 +1848,143 @@
const double Crz = orientation[8];
// Sensitivity matrix using only diagonal of A to approximate its inverse
- const double Spp = Cpx * Cpx * Aixjx + Cpy * Cpy * Aiyjy + Cpz * Cpz * Aizjz;
- const double Spq = Cpx * Cqx * Aixjx + Cpy * Cqy * Aiyjy + Cpz * Cqz * Aizjz;
- const double Spr = Cpx * Crx * Aixjx + Cpy * Cry * Aiyjy + Cpz * Crz * Aizjz;
- const double Sqq = Cqx * Cqx * Aixjx + Cqy * Cqy * Aiyjy + Cqz * Cqz * Aizjz;
- const double Sqr = Cqx * Crx * Aixjx + Cqy * Cry * Aiyjy + Cqz * Crz * Aizjz;
- const double Srr = Crx * Crx * Aixjx + Cry * Cry * Aiyjy + Crz * Crz * Aizjz;
-
- const double slipShearMag = sqrt((*slip)[0] * (*slip)[0] +
- (*slip)[1] * (*slip)[1]);
+ const double Spp = Cpx * Cpx * Aixjx + Cpy * Cpy * Aiyjy + Cpz * Cpz * Aizjz;
+ const double Spq = Cpx * Cqx * Aixjx + Cpy * Cqy * Aiyjy + Cpz * Cqz * Aizjz;
+ const double Spr = Cpx * Crx * Aixjx + Cpy * Cry * Aiyjy + Cpz * Crz * Aizjz;
+ const double Sqq = Cqx * Cqx * Aixjx + Cqy * Cqy * Aiyjy + Cqz * Cqz * Aizjz;
+ const double Sqr = Cqx * Crx * Aixjx + Cqy * Cry * Aiyjy + Cqz * Crz * Aizjz;
+ const double Srr = Crx * Crx * Aixjx + Cry * Cry * Aiyjy + Crz * Crz * Aizjz;
+
+ const double dlp = dLagrangeTpdt[0];
+ const double dlq = dLagrangeTpdt[1];
+ const double dlr = dLagrangeTpdt[2];
+ (*slip)[0] -= Spp * dlp + Spq * dlq + Spr * dlr;
+ (*slip)[1] -= Spq * dlp + Sqq * dlq + Sqr * dlr;
+ (*slip)[2] -= Spr * dlp + Sqr * dlq + Srr * dlr;
+
+ std::cout << "Slip: (" << (*slip)[0] << ", " << (*slip)[1] << ", " << (*slip)[2] << ")" << std::endl;
+
+ PetscLogFlops(3*3 + 6*8 + 3*6);
+} // _sensitivitySolveLumped3D
+
+// ----------------------------------------------------------------------
+// Constrain solution space with lumped Jacobian in 1-D.
+void
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D(double_array* dLagrangeTpdt,
+ const double_array& slip,
+ const double_array& sliprate,
+ const double_array& tractionTpdt,
+ const double area)
+{ // _constrainSolnSpace1D
+ assert(0 != dLagrangeTpdt);
+
+ if (tractionTpdt[0] < 0) {
+ // if compression, then no changes to solution
+ } else {
+ // if tension, then traction is zero.
+
+ const double dlp = -tractionTpdt[0] * area;
+ (*dLagrangeTpdt)[0] = dlp;
+ } // else
+
+ PetscLogFlops(2);
+} // _constrainSolnSpace1D
+
+// ----------------------------------------------------------------------
+// Constrain solution space with lumped Jacobian in 2-D.
+void
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D(double_array* dLagrangeTpdt,
+ const double_array& slip,
+ const double_array& slipRate,
+ const double_array& tractionTpdt,
+ const double area)
+{ // _constrainSolnSpace2D
+ assert(0 != dLagrangeTpdt);
+
+ std::cout << "Normal traction:" << tractionTpdt[1] << std::endl;
+
+ const double slipMag = fabs(slip[0]);
+ const double slipRateMag = fabs(slipRate[0]);
+
+ const double tractionNormal = tractionTpdt[1];
+ const double tractionShearMag = fabs(tractionTpdt[0]);
+
+ if (tractionNormal < 0 && 0.0 == slip[1]) {
+ // if in compression and no opening
+ std::cout << "FAULT IN COMPRESSION" << std::endl;
+ const double frictionStress = _friction->calcFriction(slipMag, slipRateMag,
+ tractionNormal);
+ std::cout << "frictionStress: " << frictionStress << std::endl;
+ if (tractionShearMag > frictionStress ||
+ (tractionShearMag < frictionStress && slipMag > 0.0)) {
+ // traction is limited by friction, so have sliding
+ std::cout << "LIMIT TRACTION, HAVE SLIDING" << std::endl;
+
+ // Update slip based on value required to stick versus friction
+ const double dlp = -(tractionShearMag - frictionStress) * area *
+ tractionTpdt[0] / tractionShearMag;
+ (*dLagrangeTpdt)[0] = dlp;
+ (*dLagrangeTpdt)[1] = 0.0;
+ } else {
+ // else friction exceeds value necessary, so stick
+ std::cout << "STICK" << std::endl;
+ // no changes to solution
+ } // if/else
+ } else {
+ // if in tension, then traction is zero.
+ std::cout << "FAULT IN TENSION" << std::endl;
+
+ (*dLagrangeTpdt)[0] = -tractionTpdt[0] * area;
+ (*dLagrangeTpdt)[1] = -tractionTpdt[1] * area;
+ } // else
+
+ PetscLogFlops(8);
+} // _constrainSolnSpace2D
+
+// ----------------------------------------------------------------------
+// Constrain solution space with lumped Jacobian in 3-D.
+void
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D(double_array* dLagrangeTpdt,
+ const double_array& slip,
+ const double_array& slipRate,
+ const double_array& tractionTpdt,
+ const double area)
+{ // _constrainSolnSpace3D
+ assert(0 != dLagrangeTpdt);
+
+ std::cout << "Normal traction:" << tractionTpdt[2] << std::endl;
+
+ const double slipShearMag = sqrt(slip[0] * slip[0] +
+ slip[1] * slip[1]);
double slipRateMag = sqrt(slipRate[0]*slipRate[0] +
- slipRate[1]*slipRate[1]);
+ slipRate[1]*slipRate[1]);
const double tractionNormal = tractionTpdt[2];
const double tractionShearMag =
sqrt(tractionTpdt[0] * tractionTpdt[0] +
- tractionTpdt[1] * tractionTpdt[1]);
+ tractionTpdt[1] * tractionTpdt[1]);
- if (tractionNormal < 0.0 && 0.0 == (*slip)[2]) {
+ if (tractionNormal < 0.0 && 0.0 == slip[2]) {
// if in compression and no opening
std::cout << "FAULT IN COMPRESSION" << std::endl;
const double frictionStress =
_friction->calcFriction(slipShearMag, slipRateMag, tractionNormal);
std::cout << "frictionStress: " << frictionStress << std::endl;
if (tractionShearMag > frictionStress ||
- (tractionShearMag < frictionStress && slipShearMag > 0.0)) {
+ (tractionShearMag < frictionStress && slipShearMag > 0.0)) {
// traction is limited by friction, so have sliding
std::cout << "LIMIT TRACTION, HAVE SLIDING" << std::endl;
// Update slip based on value required to stick versus friction
- const double dlp = (tractionShearMag - frictionStress) * area *
- tractionTpdt[0] / tractionShearMag;
- const double dlq = (tractionShearMag - frictionStress) * area *
- tractionTpdt[1] / tractionShearMag;
+ const double dlp = -(tractionShearMag - frictionStress) * area *
+ tractionTpdt[0] / tractionShearMag;
+ const double dlq = -(tractionShearMag - frictionStress) * area *
+ tractionTpdt[1] / tractionShearMag;
- (*dLagrangeTpdt)[0] = -dlp;
- (*dLagrangeTpdt)[1] = -dlq;
- (*slip)[0] += Spp * dlp + Spq * dlq;
- (*slip)[1] += Spq * dlp + Sqq * dlq;
+ (*dLagrangeTpdt)[0] = dlp;
+ (*dLagrangeTpdt)[1] = dlq;
+ (*dLagrangeTpdt)[2] = 0.0;
- std::cout << "Estimated slip: " << " " << (*slip)[0] << " "
- << (*slip)[1] << " " << (*slip)[2] << std::endl;
} else {
// else friction exceeds value necessary, so stick
std::cout << "STICK" << std::endl;
@@ -1645,383 +1994,13 @@
// if in tension, then traction is zero.
std::cout << "FAULT IN TENSION" << std::endl;
- // Update slip based on value required to stick versus
- // zero traction
- const double dlp = tractionTpdt[0] * area;
- const double dlq = tractionTpdt[1] * area;
- const double dlr = tractionTpdt[2] * area;
-
- (*dLagrangeTpdt)[0] = -dlp;
- (*dLagrangeTpdt)[1] = -dlq;
- (*dLagrangeTpdt)[2] = -dlr;
- (*slip)[0] += Spp * dlp + Spq * dlq + Spr * dlr;
- (*slip)[1] += Spq * dlp + Sqq * dlq + Sqr * dlr;
- (*slip)[2] += Spr * dlp + Sqr * dlq + Srr * dlr;
-
- std::cout << "Estimated slip: " << " " << (*slip)[0] << " "
- << (*slip)[1] << " " << (*slip)[2] << std::endl;
-
+ (*dLagrangeTpdt)[0] = -tractionTpdt[0] * area;
+ (*dLagrangeTpdt)[1] = -tractionTpdt[1] * area;
+ (*dLagrangeTpdt)[2] = -tractionTpdt[2] * area;
} // else
- PetscLogFlops(0); // :TODO: Fix this
-} // constrainSolnSpace3D
+ PetscLogFlops(22);
+} // _constrainSolnSpace3D
-// ----------------------------------------------------------------------
-// Constrain solution space in 2-D.
-void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped2x22D(
- double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& slipRate,
- const double_array& tractionTpdt,
- const double_array& orientation,
- const double_array& jacobianN,
- const double_array& jacobianP,
- const double area)
-{ // constrainSolnSpace2x22D
- assert(0 != dLagrangeTpdt);
- assert(0 != slip);
- std::cout << "Normal traction:" << tractionTpdt[1] << std::endl;
-
- // Sensitivity of slip to changes in the Lagrange multipliers
- assert(jacobianN[0] > 0.0);
- assert(jacobianP[0] > 0.0);
- assert(jacobianN[1] > 0.0);
- assert(jacobianP[1] > 0.0);
- assert(jacobianN[2] > 0.0);
- assert(jacobianP[2] > 0.0);
- assert(jacobianN[3] > 0.0);
- assert(jacobianP[3] > 0.0);
-
- // Fault orientation matrix
- const double Cpx = orientation[0];
- const double Cpy = orientation[1];
- const double Cqx = orientation[2];
- const double Cqy = orientation[3];
-
- // 3 x 3 Jacobian block of i side of fault
- const double Aixx = jacobianN[0];
- const double Aixy = jacobianN[1];
- const double Aiyx = jacobianN[2];
- const double Aiyy = jacobianN[3];
-
- // 3 x 3 Jacobian block of j side of fault
- const double Ajxx = jacobianP[0];
- const double Ajxy = jacobianP[1];
- const double Ajyx = jacobianP[2];
- const double Ajyy = jacobianP[3];
-
- // Determinant of 2 X 2 block of Jacobian on each side of the fault
- const double Deti = Aixx * Aiyy - Aixy*Aiyx;
- const double Detj = Ajxx * Ajyy - Ajxy*Ajyx;
- assert(Deti > 0.0);
- assert(Detj > 0.0);
-
- // Co-factor matrix for i side of fault
- const double Ci11 = Aiyy;
- const double Ci12 = -Aiyx;
- const double Ci21 = -Aixy;
- const double Ci22 = Aixx;
-
- // Co-factor matrix for j side of fault
- const double Cj11 = Ajyy;
- const double Cj12 = -Ajyx;
- const double Cj21 = -Ajxy;
- const double Cj22 = Ajxx;
-
- // Contribution to sensitivity from i side of fault
- const double Sipp = Ci11 * Cpx * Cpx + Ci22 * Cpy * Cpy + (Ci12 + Ci21) * Cpx * Cpy;
- const double Sipq = (Ci11 * Cpx + Ci12 * Cpy) * Cqx + (Ci21 * Cpx + Ci22 * Cpy) * Cqy;
- const double Siqp = (Ci11 * Cqx + Ci12 * Cqy) * Cpx + (Ci21 * Cqx + Ci22 * Cqy) * Cpy;
- const double Siqq = Ci11 * Cqx * Cqx + Ci22 * Cqy * Cqy + (Ci12 + Ci21) * Cqx * Cqy;
-
- // Contribution to sensitivity from i side of fault
- const double Sjpp = Cj11 * Cpx * Cpx + Cj22 * Cpy * Cpy + (Cj12 + Cj21) * Cpx * Cpy;
- const double Sjpq = (Cj11 * Cpx + Cj12 * Cpy) * Cqx + (Cj21 * Cpx + Cj22 * Cpy) * Cqy;
- const double Sjqp = (Cj11 * Cqx + Cj12 * Cqy) * Cpx + (Cj21 * Cqx + Cj22 * Cqy) * Cpy;
- const double Sjqq = Cj11 * Cqx * Cqx + Cj22 * Cqy * Cqy + (Cj12 + Cj21) * Cqx * Cqy;
-
- // Sensitivity Matrix
- const double Spp = Sipp/Deti + Sjpp/Detj;
- const double Spq = Sipq/Deti + Sjpq/Detj;
- const double Sqp = Siqp/Deti + Sjqp/Detj;
- const double Sqq = Siqq/Deti + Sjqq/Detj;
-
- const double slipMag = fabs((*slip)[0]);
- const double slipRateMag = fabs(slipRate[0]);
-
- const double tractionNormal = tractionTpdt[1];
- const double tractionShearMag = fabs(tractionTpdt[0]);
-
- if (tractionNormal < 0 && 0.0 == (*slip)[1]) {
- // if in compression and no opening
- std::cout << "FAULT IN COMPRESSION" << std::endl;
- const double frictionStress = _friction->calcFriction(slipMag, slipRateMag,
- tractionNormal);
- std::cout << "frictionStress: " << frictionStress << std::endl;
- if (tractionShearMag > frictionStress ||
- (tractionShearMag < frictionStress && slipMag > 0.0)) {
- // traction is limited by friction, so have sliding
- std::cout << "LIMIT TRACTION, HAVE SLIDING" << std::endl;
-
- // Update slip based on value required to stick versus friction
- const double dlp = (tractionShearMag - frictionStress) * area *
- tractionTpdt[0] / tractionShearMag;
- (*dLagrangeTpdt)[0] = -dlp;
- (*slip)[0] += Spp * dlp;
- std::cout << "Estimated slip: " << (*slip)[0] << std::endl;
- } else {
- // else friction exceeds value necessary, so stick
- std::cout << "STICK" << std::endl;
- // no changes to solution
- } // if/else
- } else {
- // if in tension, then traction is zero.
- std::cout << "FAULT IN TENSION" << std::endl;
-
- // Update slip based on value required to stick versus
- // zero traction
- const double dlp = tractionTpdt[0] * area;
- const double dlq = tractionTpdt[1] * area;
-
- (*dLagrangeTpdt)[0] = -dlp;
- (*dLagrangeTpdt)[1] = -dlq;
- (*slip)[0] += Spp * dlp + Spq * dlq;
- (*slip)[1] += Spq * dlp + Sqq * dlq;
- } // else
-
- PetscLogFlops(0); // :TODO: Fix this
-} // constrainSolnSpace2x22D
-
-// ----------------------------------------------------------------------
-// Constrain solution space in 3-D.
-void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceLumped3x33D(
- double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& slipRate,
- const double_array& tractionTpdt,
- const double_array& orientation,
- const double_array& jacobianN,
- const double_array& jacobianP,
- const double area)
-{ // constrainSolnSpace3x33D
- assert(0 != dLagrangeTpdt);
- assert(0 != slip);
-
- std::cout << "Normal traction:" << tractionTpdt[2] << std::endl;
-
- // Sensitivity of slip to changes in the Lagrange multipliers
- assert(jacobianN[0] > 0.0);
- assert(jacobianP[0] > 0.0);
- assert(jacobianN[1] > 0.0);
- assert(jacobianP[1] > 0.0);
- assert(jacobianN[2] > 0.0);
- assert(jacobianP[2] > 0.0);
- assert(jacobianN[3] > 0.0);
- assert(jacobianP[3] > 0.0);
- assert(jacobianN[4] > 0.0);
- assert(jacobianP[4] > 0.0);
- assert(jacobianN[5] > 0.0);
- assert(jacobianP[5] > 0.0);
- assert(jacobianN[6] > 0.0);
- assert(jacobianP[6] > 0.0);
- assert(jacobianN[7] > 0.0);
- assert(jacobianP[7] > 0.0);
- assert(jacobianN[8] > 0.0);
- assert(jacobianP[8] > 0.0);
-
- // 3 x 3 Jacobian block of i side of fault
- const double Aixx = jacobianN[0];
- const double Aixy = jacobianN[1];
- const double Aixz = jacobianN[2];
- const double Aiyx = jacobianN[3];
- const double Aiyy = jacobianN[4];
- const double Aiyz = jacobianN[5];
- const double Aizx = jacobianN[6];
- const double Aizy = jacobianN[7];
- const double Aizz = jacobianN[8];
-
- // 3 x 3 Jacobian block of j side of fault
- const double Ajxx = jacobianP[0];
- const double Ajxy = jacobianP[1];
- const double Ajxz = jacobianP[2];
- const double Ajyx = jacobianP[3];
- const double Ajyy = jacobianP[4];
- const double Ajyz = jacobianP[5];
- const double Ajzx = jacobianP[6];
- const double Ajzy = jacobianP[7];
- const double Ajzz = jacobianP[8];
-
- // Fault orientation matrix
- const double Cpx = orientation[0];
- const double Cpy = orientation[1];
- const double Cpz = orientation[2];
- const double Cqx = orientation[3];
- const double Cqy = orientation[4];
- const double Cqz = orientation[5];
- const double Crx = orientation[6];
- const double Cry = orientation[7];
- const double Crz = orientation[8];
-
- // Determinant of 3 X 3 block of Jacobian on each side of the fault
- const double Deti = Aixz * (-Aiyy * Aizx + Aiyx * Aizy) +
- Aixy * (Aiyz * Aizx - Aiyx * Aizz) + Aixx * (-Aiyz * Aizy + Aiyy * Aizz);
- const double Detj = Ajxz * (-Ajyy * Ajzx + Ajyx * Ajzy) +
- Ajxy * (Ajyz * Ajzx - Ajyx * Ajzz) + Ajxx * (-Ajyz * Ajzy + Ajyy * Ajzz);
- assert(Deti > 0.0);
- assert(Detj > 0.0);
-
- // Co-factor matrix for i side of fault
- const double Ci11 = Aiyz * Aizy + Aiyy * Aizz;
- const double Ci12 = Aiyz * Aizx - Aiyx * Aizz;
- const double Ci13 = -Aiyy * Aizx + Aiyx * Aizy;
- const double Ci21 = Aixz * Aizy - Aixy * Aizz;
- const double Ci22 = -Aixz * Aizx + Aixx * Aizz;
- const double Ci23 = Aixy * Aizx - Aixx * Aizy;
- const double Ci31 = Aixz * Aiyy + Aixy * Aiyz;
- const double Ci32 = Aixz * Aiyx - Aixx * Aiyz;
- const double Ci33 = -Aixy * Aiyx + Aixx * Aiyy;
-
- // Co-factor matrix for j side of fault
- const double Cj11 = Ajyz * Ajzy + Ajyy * Ajzz;
- const double Cj12 = Ajyz * Ajzx - Ajyx * Ajzz;
- const double Cj13 = -Ajyy * Ajzx + Ajyx * Ajzy;
- const double Cj21 = Ajxz * Ajzy - Ajxy * Ajzz;
- const double Cj22 = -Ajxz * Ajzx + Ajxx * Ajzz;
- const double Cj23 = Ajxy * Ajzx - Ajxx * Ajzy;
- const double Cj31 = Ajxz * Ajyy + Ajxy * Ajyz;
- const double Cj32 = Ajxz * Ajyx - Ajxx * Ajyz;
- const double Cj33 = -Ajxy * Ajyx + Ajxx * Ajyy;
-
- // Contribution to sensitivity from i side of fault
- const double Sipp = Ci11 * Cpx * Cpx + Ci22 * Cpy * Cpy + Ci33 * Cpz * Cpz +
- (Ci12 + Ci21) * Cpx * Cpy + (Ci13 + Ci31) * Cpx * Cpz + (Ci23 + Ci32) * Cpy * Cpz;
- const double Siqq = Ci11 * Cqx * Cqx + Ci22 * Cqy * Cqy + Ci33 * Cqz * Cqz +
- (Ci12 + Ci21) * Cqx * Cqy + (Ci13 + Ci31) * Cqx * Cqz + (Ci23 + Ci32) * Cqy * Cqz;
- const double Sirr = Ci11 * Crx * Crx + Ci22 * Cry * Cry + Ci33 * Crz * Crz +
- (Ci12 + Ci21) * Crx * Cry + (Ci13 + Ci31) * Crx * Crz + (Ci23 + Ci32) * Cry * Crz;
- const double Sipq = (Ci11 * Cpx + Ci12 * Cpy + Ci13 * Cpz) * Cqx +
- (Ci21 * Cpx + Ci22 * Cpy + Ci23 * Cpz) * Cqy +
- (Ci31 * Cpx + Ci32 * Cpy + Ci33 * Cpz) * Cqz;
- const double Sipr = (Ci11 * Cpx + Ci12 * Cpy + Ci13 * Cpz) * Crx +
- (Ci21 * Cpx + Ci22 * Cpy + Ci23 * Cpz) * Cry +
- (Ci31 * Cpx + Ci32 * Cpy + Ci33 * Cpz) * Crz;
- const double Siqp = (Ci11 * Cqx + Ci12 * Cqy + Ci13 * Cqz) * Cpx +
- (Ci21 * Cqx + Ci22 * Cqy + Ci23 * Cqz) * Cpy +
- (Ci31 * Cqx + Ci32 * Cqy + Ci33 * Cqz) * Cpz;
- const double Siqr = (Ci11 * Cqx + Ci12 * Cqy + Ci13 * Cqz) * Crx +
- (Ci21 * Cqx + Ci22 * Cqy + Ci23 * Cqz) * Cry +
- (Ci31 * Cqx + Ci32 * Cqy + Ci33 * Cqz) * Crz;
- const double Sirp = (Ci11 * Crx + Ci12 * Cry + Ci13 * Crz) * Cpx +
- (Ci21 * Crx + Ci22 * Cry + Ci23 * Crz) * Cpy +
- (Ci31 * Crx + Ci32 * Cry + Ci33 * Crz) * Cpz;
- const double Sirq = (Ci11 * Crx + Ci12 * Cry + Ci13 * Crz) * Cqx +
- (Ci21 * Crx + Ci22 * Cry + Ci23 * Crz) * Cqy +
- (Ci31 * Crx + Ci32 * Cry + Ci33 * Crz) * Cqz;
-
- // Contribution to sensitivity from i side of fault
- const double Sjpp = Cj11 * Cpx * Cpx + Cj22 * Cpy * Cpy + Cj33 * Cpz * Cpz +
- (Cj12 + Cj21) * Cpx * Cpy + (Cj13 + Cj31) * Cpx * Cpz + (Cj23 + Cj32) * Cpy * Cpz;
- const double Sjqq = Cj11 * Cqx * Cqx + Cj22 * Cqy * Cqy + Cj33 * Cqz * Cqz +
- (Cj12 + Cj21) * Cqx * Cqy + (Cj13 + Cj31) * Cqx * Cqz + (Cj23 + Cj32) * Cqy * Cqz;
- const double Sjrr = Cj11 * Crx * Crx + Cj22 * Cry * Cry + Cj33 * Crz * Crz +
- (Cj12 + Cj21) * Crx * Cry + (Cj13 + Cj31) * Crx * Crz + (Cj23 + Cj32) * Cry * Crz;
- const double Sjpq = (Cj11 * Cpx + Cj12 * Cpy + Cj13 * Cpz) * Cqx +
- (Cj21 * Cpx + Cj22 * Cpy + Cj23 * Cpz) * Cqy +
- (Cj31 * Cpx + Cj32 * Cpy + Cj33 * Cpz) * Cqz;
- const double Sjpr = (Cj11 * Cpx + Cj12 * Cpy + Cj13 * Cpz) * Crx +
- (Cj21 * Cpx + Cj22 * Cpy + Cj23 * Cpz) * Cry +
- (Cj31 * Cpx + Cj32 * Cpy + Cj33 * Cpz) * Crz;
- const double Sjqp = (Cj11 * Cqx + Cj12 * Cqy + Cj13 * Cqz) * Cpx +
- (Cj21 * Cqx + Cj22 * Cqy + Cj23 * Cqz) * Cpy +
- (Cj31 * Cqx + Cj32 * Cqy + Cj33 * Cqz) * Cpz;
- const double Sjqr = (Cj11 * Cqx + Cj12 * Cqy + Cj13 * Cqz) * Crx +
- (Cj21 * Cqx + Cj22 * Cqy + Cj23 * Cqz) * Cry +
- (Cj31 * Cqx + Cj32 * Cqy + Cj33 * Cqz) * Crz;
- const double Sjrp = (Cj11 * Crx + Cj12 * Cry + Cj13 * Crz) * Cpx +
- (Cj21 * Crx + Cj22 * Cry + Cj23 * Crz) * Cpy +
- (Cj31 * Crx + Cj32 * Cry + Cj33 * Crz) * Cpz;
- const double Sjrq = (Cj11 * Crx + Cj12 * Cry + Cj13 * Crz) * Cqx +
- (Cj21 * Crx + Cj22 * Cry + Cj23 * Crz) * Cqy +
- (Cj31 * Crx + Cj32 * Cry + Cj33 * Crz) * Cqz;
-
- // Sensitivity Matrix
- const double Spp = Sipp/Deti + Sjpp/Detj;
- const double Spq = Sipq/Deti + Sjpq/Detj;
- const double Spr = Sipr/Deti + Sjpr/Detj;
- const double Sqp = Siqp/Deti + Sjqp/Detj;
- const double Sqq = Siqq/Deti + Sjqq/Detj;
- const double Sqr = Siqr/Deti + Sjqr/Detj;
- const double Srp = Sirp/Deti + Sjrp/Detj;
- const double Srq = Sirq/Deti + Sjrq/Detj;
- const double Srr = Sirr/Deti + Sjrr/Detj;
-
-
- const double slipShearMag = sqrt((*slip)[0] * (*slip)[0] +
- (*slip)[1] * (*slip)[1]);
- double slipRateMag = sqrt(slipRate[0]*slipRate[0] +
- slipRate[1]*slipRate[1]);
-
- const double tractionNormal = tractionTpdt[2];
- const double tractionShearMag =
- sqrt(tractionTpdt[0] * tractionTpdt[0] +
- tractionTpdt[1] * tractionTpdt[1]);
-
- if (tractionNormal < 0.0 && 0.0 == (*slip)[2]) {
- // if in compression and no opening
- std::cout << "FAULT IN COMPRESSION" << std::endl;
- const double frictionStress =
- _friction->calcFriction(slipShearMag, slipRateMag, tractionNormal);
- std::cout << "frictionStress: " << frictionStress << std::endl;
- if (tractionShearMag > frictionStress ||
- (tractionShearMag < frictionStress && slipShearMag > 0.0)) {
- // traction is limited by friction, so have sliding
- std::cout << "LIMIT TRACTION, HAVE SLIDING" << std::endl;
-
- // Update slip based on value required to stick versus friction
- const double dlp = (tractionShearMag - frictionStress) * area *
- tractionTpdt[0] / tractionShearMag;
- const double dlq = (tractionShearMag - frictionStress) * area *
- tractionTpdt[1] / tractionShearMag;
-
- (*dLagrangeTpdt)[0] = -dlp;
- (*dLagrangeTpdt)[1] = -dlq;
- (*slip)[0] += Spp * dlp + Spq * dlq;
- (*slip)[1] += Spq * dlp + Sqq * dlq;
-
- std::cout << "Estimated slip: " << " " << (*slip)[0] << " "
- << (*slip)[1] << " " << (*slip)[2] << std::endl;
- } else {
- // else friction exceeds value necessary, so stick
- std::cout << "STICK" << std::endl;
- // no changes to solution
- } // if/else
- } else {
- // if in tension, then traction is zero.
- std::cout << "FAULT IN TENSION" << std::endl;
-
- // Update slip based on value required to stick versus
- // zero traction
- const double dlp = tractionTpdt[0] * area;
- const double dlq = tractionTpdt[1] * area;
- const double dlr = tractionTpdt[2] * area;
-
- (*dLagrangeTpdt)[0] = -dlp;
- (*dLagrangeTpdt)[1] = -dlq;
- (*dLagrangeTpdt)[2] = -dlr;
- (*slip)[0] += Spp * dlp + Spq * dlq + Spr * dlr;
- (*slip)[1] += Spq * dlp + Sqq * dlq + Sqr * dlr;
- (*slip)[2] += Spr * dlp + Sqr * dlq + Srr * dlr;
-
- std::cout << "Estimated slip: " << " " << (*slip)[0] << " "
- << (*slip)[1] << " " << (*slip)[2] << std::endl;
-
- } // else
-
- PetscLogFlops(0); // :TODO: Fix this
-} // constrainSolnSpace3x33D
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2010-03-23 22:52:56 UTC (rev 16443)
@@ -23,6 +23,7 @@
#include "FaultCohesiveLagrange.hh" // ISA FaultCohesiveLagrange
#include "pylith/friction/frictionfwd.hh" // HOLDSA Friction model
+#include "pylith/utils/petscfwd.h" // HASA PetscKSP
// FaultCohesiveDyn -----------------------------------------------------
/**
@@ -166,115 +167,128 @@
*/
void _updateSlipRate(const topology::SolutionFields& fields);
- /** Update Jacobian blocks associated with DOF of vertices on
- * negative and positive sides of the fault associated with Lagrange
- * vertex k.
+ /** Setup sensitivity problem to compute change in slip given change
+ * in Lagrange multipliers.
*
- * @param jacobian Jacobian matrix.
+ * @param jacobian Jacobian matrix for entire domain.
+ */
+ void _sensitivitySetup(const topology::Jacobian& jacobian);
+
+ /** Update the Jacobian values for the sensitivity solve.
+ *
+ * @param negativeSide True if solving sensitivity problem for
+ * negative side of the fault, false if solving sensitivity problem
+ * for positive side of the fault.
+ * @param jacobian Jacobian matrix for entire domain.
* @param fields Solution fields.
*/
- void _updateJacobianBlocks(const topology::Jacobian& jacobian,
- const topology::SolutionFields& fields);
+ void _sensitivityUpdateJacobian(const bool negativeSide,
+ const topology::Jacobian& jacobian,
+ const topology::SolutionFields& fields);
- /** Constrain solution space with lumped Jacobian in 1-D.
+ /** Reform residual for sensitivity problem.
*
- * @param dLagrangeTpdt Adjustment to Lagrange multiplier.
+ * @param negativeSide True if solving sensitivity problem for
+ * negative side of the fault, false if solving sensitivity problem
+ * for positive side of the fault.
+ */
+ void _sensitivityReformResidual(const bool negativeSide);
+
+ /// Solve sensitivity problem.
+ void _sensitivitySolve(void);
+
+ /** Update the solution (displacement increment) values based on
+ * the sensitivity solve.
+ *
+ * @param negativeSide True if solving sensitivity problem for
+ * negative side of the fault, false if solving sensitivity problem
+ * for positive side of the fault.
+ */
+ void _sensitivityUpdateSoln(const bool negativeSide);
+
+ /** Solve slip/Lagrange multiplier sensitivity problem for case of lumped Jacobian in 1-D.
+ *
* @param slip Adjustment to slip assoc. w/Lagrange multiplier vertex.
- * @param slipRate Slip rate assoc. w/Lagrange multiplier vertex.
- * @param tractionTpdt Fault traction assoc. w/Lagrange multiplier vertex.
+ * @param dLagrangeTpdt Change in Lagrange multiplier.
* @param orientation Fault orientation assoc. w/Lagrang multiplier vertex.
* @param jacobianN Jacobian for vertex on - side of the fault.
* @param jacobianP Jacobian for vertex on + side of the fault.
- * @param area Fault area associated w/Lagrange multiplier vertex.
*/
- void _constrainSolnSpaceLumped1D(double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& sliprate,
- const double_array& tractionTpdt,
- const double_array& orientation,
- const double_array& jacobianN,
- const double_array& jacobianP,
- const double area);
+ void _sensitivitySolveLumped1D(double_array* slip,
+ const double_array& dLagrangeTpdt,
+ const double_array& orientation,
+ const double_array& jacobianN,
+ const double_array& jacobianP);
- /** Constrain solution space with lumped Jacobian in 2-D.
+ /** Solve slip/Lagrange multiplier sensitivity problem for case of lumped Jacobian in 2-D.
*
- * @param dLagrangeTpdt Adjustment to Lagrange multiplier.
* @param slip Adjustment to slip assoc. w/Lagrange multiplier vertex.
- * @param slipRate Slip rate assoc. w/Lagrange multiplier vertex.
- * @param tractionTpdt Fault traction assoc. w/Lagrange multiplier vertex.
+ * @param dLagrangeTpdt Change in Lagrange multiplier.
* @param orientation Fault orientation assoc. w/Lagrang multiplier vertex.
* @param jacobianN Jacobian for vertex on - side of the fault.
* @param jacobianP Jacobian for vertex on + side of the fault.
- * @param area Fault area associated w/Lagrange multiplier vertex.
*/
- void _constrainSolnSpaceLumped2D(double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& sliprate,
- const double_array& tractionTpdt,
- const double_array& orientation,
- const double_array& jacobianN,
- const double_array& jacobianP,
- const double area);
+ void _sensitivitySolveLumped2D(double_array* slip,
+ const double_array& dLagrangeTpdt,
+ const double_array& orientation,
+ const double_array& jacobianN,
+ const double_array& jacobianP);
- /** Constrain solution space with lumped Jacobian in 3-D.
+ /** Solve slip/Lagrange multiplier sensitivity problem for case of lumped Jacobian in 3-D.
*
- * @param dLagrangeTpdt Adjustment to Lagrange multiplier.
* @param slip Adjustment to slip assoc. w/Lagrange multiplier vertex.
- * @param slipRate Slip rate assoc. w/Lagrange multiplier vertex.
- * @param tractionTpdt Fault traction assoc. w/Lagrange multiplier vertex.
+ * @param dLagrangeTpdt Change in Lagrange multiplier.
* @param orientation Fault orientation assoc. w/Lagrang multiplier vertex.
* @param jacobianN Jacobian for vertex on - side of the fault.
* @param jacobianP Jacobian for vertex on + side of the fault.
+ */
+ void _sensitivitySolveLumped3D(double_array* slip,
+ const double_array& dLagrangeTpdt,
+ const double_array& orientation,
+ const double_array& jacobianN,
+ const double_array& jacobianP);
+
+ /** Constrain solution space with lumped Jacobian in 1-D.
+ *
+ * @param dLagrangeTpdt Adjustment to Lagrange multiplier.
+ * @param slip Slip assoc. w/Lagrange multiplier vertex.
+ * @param slipRate Slip rate assoc. w/Lagrange multiplier vertex.
+ * @param tractionTpdt Fault traction assoc. w/Lagrange multiplier vertex.
* @param area Fault area associated w/Lagrange multiplier vertex.
*/
- void _constrainSolnSpaceLumped3D(double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& sliprate,
- const double_array& tractionTpdt,
- const double_array& orientation,
- const double_array& jacobianN,
- const double_array& jacobianP,
- const double area);
+ void _constrainSolnSpace1D(double_array* dLagrangeTpdt,
+ const double_array& slip,
+ const double_array& slipRate,
+ const double_array& tractionTpdt,
+ const double area);
- /** Constrain solution space with 2 x 2 Jacobian block in 2-D.
+ /** Constrain solution space with lumped Jacobian in 2-D.
*
* @param dLagrangeTpdt Adjustment to Lagrange multiplier.
- * @param slip Adjustment to slip assoc. w/Lagrange multiplier vertex.
+ * @param slip Slip assoc. w/Lagrange multiplier vertex.
* @param slipRate Slip rate assoc. w/Lagrange multiplier vertex.
* @param tractionTpdt Fault traction assoc. w/Lagrange multiplier vertex.
- * @param orientation Fault orientation assoc. w/Lagrang multiplier vertex.
- * @param jacobianN Jacobian for vertex on - side of the fault.
- * @param jacobianP Jacobian for vertex on + side of the fault.
* @param area Fault area associated w/Lagrange multiplier vertex.
*/
- void _constrainSolnSpaceLumped2x22D(double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& sliprate,
- const double_array& tractionTpdt,
- const double_array& orientation,
- const double_array& jacobianN,
- const double_array& jacobianP,
- const double area);
+ void _constrainSolnSpace2D(double_array* dLagrangeTpdt,
+ const double_array& slip,
+ const double_array& slipRate,
+ const double_array& tractionTpdt,
+ const double area);
- /** Constrain solution space with 3 X 3 Jacobian block in 3-D.
+ /** Constrain solution space with lumped Jacobian in 3-D.
*
* @param dLagrangeTpdt Adjustment to Lagrange multiplier.
- * @param slip Adjustment to slip assoc. w/Lagrange multiplier vertex.
+ * @param slip Slip assoc. w/Lagrange multiplier vertex.
* @param slipRate Slip rate assoc. w/Lagrange multiplier vertex.
* @param tractionTpdt Fault traction assoc. w/Lagrange multiplier vertex.
- * @param orientation Fault orientation assoc. w/Lagrang multiplier vertex.
- * @param jacobianN Jacobian for vertex on - side of the fault.
- * @param jacobianP Jacobian for vertex on + side of the fault.
* @param area Fault area associated w/Lagrange multiplier vertex.
*/
- void _constrainSolnSpaceLumped3x33D(double_array* dLagrangeTpdt,
- double_array* slip,
- const double_array& sliprate,
- const double_array& tractionTpdt,
- const double_array& orientation,
- const double_array& jacobianN,
- const double_array& jacobianP,
- const double area);
+ void _constrainSolnSpace3D(double_array* dLagrangeTpdt,
+ const double_array& slip,
+ const double_array& slipRate,
+ const double_array& tractionTpdt,
+ const double area);
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
@@ -285,7 +299,12 @@
/// To identify constitutive model
friction::FrictionModel* _friction;
- // NOT IMPLEMENTED ////////////////////////////////////////////////////
+ /// Sparse matrix for sensitivity solve.
+ topology::Jacobian* _jacobian;
+
+ PetscKSP _ksp; ///< PETSc KSP linear solver for sensitivity problem.
+
+// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
/// Not implemented
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-03-23 22:52:56 UTC (rev 16443)
@@ -87,7 +87,7 @@
delete _faultMesh;
_faultMesh = new topology::SubMesh();
- CohesiveTopology::createFaultParallel(_faultMesh, mesh, id(),
+ CohesiveTopology::createFaultParallel(_faultMesh, mesh, id(),
_useLagrangeConstraints);
_initializeCohesiveInfo(mesh);
@@ -763,12 +763,18 @@
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& faultCells =
+ faultSieveMesh->heightStratum(0);
+ assert(!faultCells.isNull());
+ SieveSubMesh::label_sequence::iterator f_iter = faultCells->begin();
+
SubMesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
const SieveSubMesh::renumbering_type::const_iterator renumberingEnd =
renumbering.end();
const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
faultSieveMesh->depthStratum(0);
+ _cohesiveToFault.clear();
typedef std::map<int, int> indexmap_type;
indexmap_type indexMap;
_cohesiveVertices.resize(vertices->size());
@@ -783,7 +789,9 @@
for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
- ++c_iter) {
+ ++c_iter, ++f_iter) {
+ _cohesiveToFault[*c_iter] = *f_iter;
+
// Get oriented closure
ncV.clear();
ALE::ISieveTraversal<SieveMesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh 2010-03-23 22:52:56 UTC (rev 16443)
@@ -291,6 +291,10 @@
/// Array of cohesive vertex information.
std::vector<CohesiveInfo> _cohesiveVertices;
+ /// Map label of cohesive cell to label of cells in fault mesh.
+ std::map<topology::Mesh::SieveMesh::point_type,
+ topology::SubMesh::SieveMesh::point_type> _cohesiveToFault;
+
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2010-03-23 22:52:56 UTC (rev 16443)
@@ -31,7 +31,6 @@
_gravityField(0),
_logger(0),
_needNewJacobian(true),
- _needJacobianDiag(false),
_needVelocity(false),
_useSolnIncr(false)
{ // constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-03-23 22:52:56 UTC (rev 16443)
@@ -108,12 +108,6 @@
virtual
bool needNewJacobian(void) const;
- /** Check whether integrator needs Jacobian diagonal.
- *
- * @returns True if integrator needs Jacobian diagonal for computation.
- */
- bool needJacobianDiag(void) const;
-
/** Check whether integrator needs velocity.
*
* @returns True if integrator needs velocity for computation.
@@ -317,10 +311,6 @@
/// Default is false;
bool _needNewJacobian;
- /// True if we need to store Jacobian diagonal, false otherwise.
- /// Default is false;
- bool _needJacobianDiag;
-
/// True if we need to compute velocity field, false otherwise.
/// Default is false;
bool _needVelocity;
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-03-23 22:52:56 UTC (rev 16443)
@@ -30,14 +30,6 @@
return _needNewJacobian;
} // needNewJacobian
-// Check whether integrator needs Jacobian diagonal.
-template<typename quadrature_type>
-inline
-bool
-pylith::feassemble::Integrator<quadrature_type>::needJacobianDiag(void) const {
- return _needJacobianDiag;
-} // needsJacobianDiag
-
// Check whether integrator needs velocity.
template<typename quadrature_type>
inline
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-03-23 22:52:56 UTC (rev 16443)
@@ -32,7 +32,6 @@
_jacobian(0),
_jacobianLumped(0),
_fields(0),
- _needJacobianDiag(false),
_needVelocity(false)
{ // constructor
} // constructor
@@ -102,24 +101,9 @@
pylith::problems::Formulation::initialize(void)
{ // initialize
- // Determine whether we need to store Jacobian diagonal
- _needJacobianDiag = false;
- int numIntegrators = _meshIntegrators.size();
- for (int i=0; i < numIntegrators; ++i) {
- assert(0 != _meshIntegrators[i]);
- if (_meshIntegrators[i]->needJacobianDiag())
- _needJacobianDiag = true;
- } // for
- numIntegrators = _submeshIntegrators.size();
- for (int i=0; i < numIntegrators; ++i) {
- assert(0 != _submeshIntegrators[i]);
- if (_submeshIntegrators[i]->needJacobianDiag())
- _needJacobianDiag = true;
- } // for
-
// Determine whether we need to compute velocity
_needVelocity = false;
- numIntegrators = _meshIntegrators.size();
+ int numIntegrators = _meshIntegrators.size();
for (int i=0; i < numIntegrators; ++i) {
assert(0 != _meshIntegrators[i]);
if (_meshIntegrators[i]->needVelocity())
@@ -324,21 +308,6 @@
// Assemble jacobian.
_jacobian->assemble("final_assembly");
- if (_needJacobianDiag) {
- if (!_fields->hasField("Jacobian diagonal")) {
- _fields->add("Jacobian diagonal", "jacobian_diagonal");
- topology::Field<topology::Mesh>& jacobianDiag = _fields->get(
- "Jacobian diagonal");
- const topology::Field<topology::Mesh>& disp = _fields->get("disp(t)");
- jacobianDiag.cloneSection(disp);
- jacobianDiag.createVector();
- jacobianDiag.createScatter();
- } // if
- topology::Field<topology::Mesh>& jacobianDiag = _fields->get(
- "Jacobian diagonal");
- MatGetDiagonal(_jacobian->matrix(), jacobianDiag.vector());
- jacobianDiag.scatterVectorToSection();
- } // if
} // reformJacobian
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-03-23 22:52:56 UTC (rev 16443)
@@ -172,7 +172,6 @@
///< Integrators over lower-dimensional subdomains of the mesh.
std::vector<IntegratorSubMesh*> _submeshIntegrators;
- bool _needJacobianDiag; ///< Integrator(s) need Jacobian diagonal.
bool _needVelocity; ///< Integrator(s) need velocity.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh 2010-03-23 22:52:56 UTC (rev 16443)
@@ -46,6 +46,7 @@
typedef Mesh::IntSection IntSection;
typedef Mesh::RestrictVisitor RestrictVisitor;
typedef Mesh::UpdateAddVisitor UpdateAddVisitor;
+ typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveMesh::order_type,PetscInt> IndicesVisitor;
// Sieve mesh for higher level domain (mesh, not submesh)
typedef Mesh::SieveMesh DomainSieveMesh;
Modified: short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/pylithapp.cfg 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/pylithapp.cfg 2010-03-23 22:52:56 UTC (rev 16443)
@@ -74,7 +74,6 @@
# ----------------------------------------------------------------------
[pylithapp.petsc]
ksp_type = gmres
-
pc_type = asm
# Change the preconditioner settings.
@@ -84,9 +83,8 @@
ksp_atol = 1.0e-15
ksp_max_it = 200
ksp_gmres_restart = 250
+snes_max_it = 300
-snes_max_it = 500
-
ksp_monitor = true
ksp_view = true
ksp_converged_reason = true
@@ -95,7 +93,7 @@
snes_view = true
snes_converged_reason = true
-log_summary = true
+#log_summary = true
# ----------------------------------------------------------------------
# output
Modified: short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/shear-sliding.cfg
===================================================================
--- short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/shear-sliding.cfg 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/shear-sliding.cfg 2010-03-23 22:52:56 UTC (rev 16443)
@@ -27,7 +27,7 @@
db_initial = spatialdata.spatialdb.UniformDB
db_initial.label = Dirichlet BC +x edge
db_initial.values = [displacement-x,displacement-y]
-db_initial.data = [0.0*m,30.0*km]
+db_initial.data = [0.0*m,32.0*km]
[pylithapp.timedependent.bc.x_neg]
Modified: short/3D/PyLith/trunk/playpen/friction/twohex8/opening.cfg
===================================================================
--- short/3D/PyLith/trunk/playpen/friction/twohex8/opening.cfg 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/playpen/friction/twohex8/opening.cfg 2010-03-23 22:52:56 UTC (rev 16443)
@@ -136,7 +136,7 @@
ksp_rtol = 1.0e-8
ksp_max_it = 100
ksp_gmres_restart = 50
-snes_max_it = 200
+snes_max_it = 10
ksp_monitor = true
ksp_view = true
@@ -145,7 +145,6 @@
snes_monitor = true
snes_view = true
snes_converged_reason = true
-snes_max_it = 1000
#log_summary = true
#start_in_debugger = true
Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py 2010-03-23 22:52:56 UTC (rev 16443)
@@ -93,7 +93,7 @@
self._info.log("Creating Jacobian matrix.")
from pylith.topology.Jacobian import Jacobian
- self.jacobian = Jacobian(self.fields)
+ self.jacobian = Jacobian(self.fields.solution(), self.matrixType)
self.jacobian.zero() # TEMPORARY, to get correct memory usage
self._debug.log(resourceUsageString())
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2010-03-23 22:52:56 UTC (rev 16443)
@@ -476,7 +476,7 @@
self._eventLogger.stagePop()
if self.viewJacobian:
- self.jacobianViewer.write(self.jacobian, t)
+ self.jacobianViewer.write(self.jacobian, t, self.mesh.comm())
self._debug.log(resourceUsageString())
return
Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py 2010-03-23 22:52:56 UTC (rev 16443)
@@ -132,7 +132,7 @@
# Allocates memory for nonzero pattern and Jacobian
self._info.log("Creating Jacobian matrix.")
from pylith.topology.Jacobian import Jacobian
- self.jacobian = Jacobian(self.fields, self.matrixType)
+ self.jacobian = Jacobian(self.fields.solution(), self.matrixType)
self.jacobian.zero() # TEMPORARY, to get correct memory usage
self._debug.log(resourceUsageString())
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2010-03-23 22:52:56 UTC (rev 16443)
@@ -221,7 +221,7 @@
const ALE::Obj<RealSection>& solutionSection = solution.section();
CPPUNIT_ASSERT(!solutionSection.isNull());
- topology::Jacobian jacobian(fields);
+ topology::Jacobian jacobian(solution);
const double t = 1.0;
bc.integrateJacobian(&jacobian, t, &fields);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2010-03-23 22:52:56 UTC (rev 16443)
@@ -77,16 +77,6 @@
} // testConstructor
// ----------------------------------------------------------------------
-// Test needJacobianDiag().
-void
-pylith::faults::TestFaultCohesiveDyn::testNeedJacobianDiag(void)
-{ // testNeedJacobianDiag
- FaultCohesiveDyn fault;
-
- CPPUNIT_ASSERT_EQUAL(true, fault.needJacobianDiag());
-} // testNeedJacobianDiag
-
-// ----------------------------------------------------------------------
// Test needVelocity().
void
pylith::faults::TestFaultCohesiveDyn::testNeedVelocity(void)
@@ -189,7 +179,7 @@
// Initial forces/tractions
if (0 != fault._dbInitialTract) {
- fault._fields->get("initial forces").view("INITIAL FORCES"); // DEBUGGING
+ //fault._fields->get("initial forces").view("INITIAL FORCES"); // DEBUGGING
const ALE::Obj<RealSection>& forcesInitialSection =
fault._fields->get("initial forces").section();
CPPUNIT_ASSERT(!forcesInitialSection.isNull());
@@ -225,7 +215,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields, "seqdense");
+ topology::Jacobian jacobian(fields.solution(), "seqdense");
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrStick);
const int spaceDim = _data->spaceDim;
@@ -331,7 +321,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields, "seqdense");
+ topology::Jacobian jacobian(fields.solution(), "seqdense");
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrSlip);
const int spaceDim = _data->spaceDim;
@@ -454,7 +444,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields, "seqdense");
+ topology::Jacobian jacobian(fields.solution(), "seqdense");
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrOpen);
const int spaceDim = _data->spaceDim;
@@ -577,7 +567,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields, "seqdense");
+ topology::Jacobian jacobian(fields.solution(), "seqdense");
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrSlip);
const int spaceDim = _data->spaceDim;
@@ -602,7 +592,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields, "seqdense");
+ topology::Jacobian jacobian(fields.solution(), "seqdense");
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrStick);
CPPUNIT_ASSERT(0 != fault._faultMesh);
@@ -835,20 +825,10 @@
rows[iRow] = iRow;
for (int iCol=0; iCol < ncols; ++iCol)
cols[iCol] = iCol;
- MatSetValues(jacobianMat, nrows, &rows[0], ncols, &cols[0], _data->jacobian, INSERT_VALUES);
+ MatSetValues(jacobianMat, nrows, &rows[0], ncols, &cols[0],
+ _data->jacobian, INSERT_VALUES);
jacobian->assemble("final_assembly");
- // Set Jacobian diagonal
- fields->add("Jacobian diagonal", "jacobian_diagonal");
- topology::Field<topology::Mesh>& jacobianDiag =
- fields->get("Jacobian diagonal");
- const topology::Field<topology::Mesh>& disp =
- fields->get("disp(t)");
- jacobianDiag.cloneSection(disp);
- jacobianDiag.createVector();
- jacobianDiag.createScatter();
- MatGetDiagonal(jacobian->matrix(), jacobianDiag.vector());
- jacobianDiag.scatterVectorToSection();
} // _setFieldsJacobian
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2010-03-23 22:52:56 UTC (rev 16443)
@@ -46,7 +46,6 @@
CPPUNIT_TEST_SUITE( TestFaultCohesiveDyn );
CPPUNIT_TEST( testConstructor );
- CPPUNIT_TEST( testNeedJacobianDiag );
CPPUNIT_TEST( testNeedVelocity );
CPPUNIT_TEST( testDBInitialTract );
@@ -82,9 +81,6 @@
/// Test constructor.
void testConstructor(void);
- /// Test needJacobianDiag().
- void testNeedJacobianDiag(void);
-
/// Test needVelocity().
void testNeedVelocity(void);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2010-03-23 22:50:46 UTC (rev 16442)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2010-03-23 22:52:56 UTC (rev 16443)
@@ -417,7 +417,7 @@
dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
} // for
- topology::Jacobian jacobian(fields);
+ topology::Jacobian jacobian(fields.solution());
const double t = 2.134;
fault.integrateJacobian(&jacobian, t, &fields);
@@ -496,7 +496,7 @@
++v_iter, ++iVertex)
dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
- topology::Jacobian jacobian(fields);
+ topology::Jacobian jacobian(fields.solution());
const double t = 2.134;
fault.integrateJacobianAssembled(&jacobian, t, &fields);
More information about the CIG-COMMITS
mailing list