[cig-commits] r18998 - in short/3D/PyLith/branches/v1.6-revisedfault: libsrc/pylith/faults unittests/libtests/faults unittests/libtests/faults/data
brad at geodynamics.org
brad at geodynamics.org
Fri Sep 30 16:02:11 PDT 2011
Author: brad
Date: 2011-09-30 16:02:10 -0700 (Fri, 30 Sep 2011)
New Revision: 18998
Modified:
short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.hh
short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveKin.cc
short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.hh
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynData.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynData.hh
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.hh
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.hh
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.hh
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.hh
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.hh
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/hex8_initialtract.spatialdb
Log:
More work on updating spontaneous rupture implementation.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -136,12 +136,13 @@
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
assert(0 != cs);
- // Create field for slip rate associated with Lagrange vertex k
- _fields->add("slip rate", "slip_rate");
- topology::Field<topology::SubMesh>& slipRate = _fields->get("slip rate");
- topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- slipRate.cloneSection(slip);
- slipRate.vectorFieldType(topology::FieldBase::VECTOR);
+ // Create field for relative velocity associated with Lagrange vertex k
+ _fields->add("relative velocity", "relative_velocity");
+ topology::Field<topology::SubMesh>& velRel =
+ _fields->get("relative velocity");
+ topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
+ velRel.cloneSection(dispRel);
+ velRel.vectorFieldType(topology::FieldBase::VECTOR);
//logger.stagePop();
} // initialize
@@ -156,6 +157,7 @@
{ // integrateResidual
assert(0 != fields);
assert(0 != _fields);
+ assert(0 != _logger);
// Initial fault tractions have been assembled, so they do not need
// assembling across processors.
@@ -166,65 +168,214 @@
if (0 == _dbInitialTract)
return;
+ const int setupEvent = _logger->eventId("FaIR setup");
+ const int geometryEvent = _logger->eventId("FaIR geometry");
+ const int computeEvent = _logger->eventId("FaIR compute");
+ const int restrictEvent = _logger->eventId("FaIR restrict");
+ const int updateEvent = _logger->eventId("FaIR update");
+
+ _logger->eventBegin(setupEvent);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ assert(cellDim == spaceDim-1);
- // Get sections
- double_array forcesInitialVertex(spaceDim);
- const ALE::Obj<RealSection>& forcesInitialSection =
- _fields->get("initial forces").section();
- assert(!forcesInitialSection.isNull());
+ // Get cohesive cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", id());
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ const int numCells = cells->size();
- double_array residualVertex(spaceDim);
+ // Get sections associated with cohesive cells
+ double_array residualCell(3*numBasis*spaceDim);
const ALE::Obj<RealSection>& residualSection = residual.section();
assert(!residualSection.isNull());
+ UpdateAddVisitor residualVisitor(*residualSection, &residualCell[0]);
- const ALE::Obj<RealSection>& slipSection = _fields->get("slip").section();
- assert(!slipSection.isNull());
+ // Get fault cell information
+ const ALE::Obj<SieveMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& faultCells =
+ faultSieveMesh->heightStratum(0);
+ assert(!faultCells.isNull());
+ assert(faultCells->size() == cells->size());
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
- residualSection);
- assert(!globalOrder.isNull());
+ // Get sections associated with fault cells
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
- const int numVertices = _cohesiveVertices.size();
- 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;
+ double_array dispRelCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispRelSection =
+ _fields->get("relative disp").section();
+ assert(!dispRelSection.isNull());
+ RestrictVisitor dispRelVisitor(*dispRelSection,
+ dispRelCell.size(), &dispRelCell[0]);
- // Compute contribution only if Lagrange constraint is local.
- if (!globalOrder->isLocal(v_lagrange))
- continue;
+ double_array initialTractionsCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& initialTractionsSection =
+ _fields->get("initial tractions").section();
+ assert(!initialTractionsSection.isNull());
+ RestrictVisitor initialTractionsVisitor(*initialTractionsSection,
+ initialTractionsCell.size(),
+ &initialTractionsCell[0]);
- // Get initial forces at fault vertex. Forces are in the global
- // coordinate system so no rotation is necessary.
- forcesInitialSection->restrictPoint(v_fault,
- &forcesInitialVertex[0],
- forcesInitialVertex.size());
+ double_array orientationCell(numBasis*spaceDim*spaceDim);
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
+ assert(!orientationSection.isNull());
+ RestrictVisitor orientationVisitor(*orientationSection,
+ orientationCell.size(),
+ &orientationCell[0]);
- assert(spaceDim == slipSection->getFiberDimension(v_fault));
- const double* slipVertex = slipSection->restrictPoint(v_fault);
- assert(0 != slipVertex);
-
- // only apply initial tractions if there is no opening
- if (0.0 == slipVertex[spaceDim-1]) {
- residualVertex = forcesInitialVertex;
+ _logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
- assert(residualVertex.size() ==
- residualSection->getFiberDimension(v_positive));
- residualSection->updateAddPoint(v_positive, &residualVertex[0]);
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ topology::SubMesh::SieveMesh::point_type c_fault =
+ _cohesiveToFault[*c_iter];
+
+ // Compute geometry information for current cell
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(geometryEvent);
+#endif
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(c_fault);
+#else
+ coordsVisitor.clear();
+ faultSieveMesh->restrictClosure(c_fault, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, c_fault);
+#endif
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(geometryEvent);
+ _logger->eventBegin(restrictEvent);
+#endif
+
+ // Restrict input fields to cell
+ initialTractionsVisitor.clear();
+ sieveMesh->restrictClosure(c_fault, initialTractionsVisitor);
+
+ dispRelVisitor.clear();
+ faultSieveMesh->restrictClosure(c_fault, dispRelVisitor);
+
+ orientationVisitor.clear();
+ faultSieveMesh->restrictClosure(c_fault, orientationVisitor);
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(restrictEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ residualCell = 0.0;
+
+ // Compute action for positive side of fault and Lagrange constraints
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const int iQ = iQuad * numBasis;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const double valI = wt*basis[iQ+iBasis];
+
+ // Add entries to residual at
+ // iN = DOF at negative node
+ // iP = DOF at positive node
+ // iL = DOF at constraint node
+
+ // Indices for negative vertex
+ const int iBN = 0*numBasis*spaceDim + iBasis*spaceDim;
+
+ // Indices for positive vertex
+ const int iBP = 1*numBasis*spaceDim + iBasis*spaceDim;
+
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double valIJ = valI * basis[iQ+jBasis];
+
+ // Indices for fault vertex
+ const int jB = jBasis*spaceDim;
+
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ // negative side of the fault
+ residualCell[iBN + iDim] -= valIJ * initialTractionsCell[jB + iDim];
+
+ // positive side of the fault
+ residualCell[iBP + iDim] += valIJ * initialTractionsCell[jB + iDim];
+
+#if 0
+ std::cout << "iBasis: " << iBasis
+ << ", jBasis: " << jBasis
+ << ", iDim: " << iDim
+ << ", valIJ: " << valIJ
+ << ", jacobianDet: " << jacobianDet[iQuad]
+ << ", residualN: " << residualCell[iBN + iDim]
+ << ", residualP: " << residualCell[iBP + iDim]
+ << std::endl;
+#endif
+
+ } // for
+ } // for
+ } // for
+ } // for
+
+
+ // Only apply initial tractions if there is no opening.
+ // If there is opening, zero out initial tractions
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const int iB = iBasis*spaceDim;
+ const int iO = iBasis*spaceDim*spaceDim;
- residualVertex *= -1.0;
- assert(residualVertex.size() ==
- residualSection->getFiberDimension(v_negative));
- residualSection->updateAddPoint(v_negative, &residualVertex[0]);
- } // if
+ double slipNormal = 0.0;
+ const int indexN = spaceDim - 1;
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
+ slipNormal += orientationCell[iO + jDim*spaceDim]*dispRelCell[iB+jDim];
+ } // for
+
+ if (slipNormal > _zeroTolerance) {
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ residualCell[iB+iDim] = 0.0;
+ } // for
+ } // if
+ } // for
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
+
+ // Assemble cell contribution into field
+ residualVisitor.clear();
+ sieveMesh->updateClosure(*c_iter, residualVisitor);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(updateEvent);
+#endif
} // for
+ PetscLogFlops(numCells*spaceDim*spaceDim*8);
- PetscLogFlops(numVertices*spaceDim);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
} // integrateResidual
// ----------------------------------------------------------------------
@@ -237,7 +388,7 @@
assert(0 != fields);
assert(0 != _fields);
- _updateSlipRate(*fields);
+ _updateVelRel(*fields);
const int spaceDim = _quadrature->spaceDim();
@@ -245,12 +396,15 @@
double_array tractionTpdtVertex(spaceDim); // Fault coordinate system
// Get sections
- const ALE::Obj<RealSection>& slipSection = _fields->get("slip").section();
- assert(!slipSection.isNull());
+ double_array slipVertex(spaceDim);
+ const ALE::Obj<RealSection>& dispRelSection =
+ _fields->get("relative disp").section();
+ assert(!dispRelSection.isNull());
- const ALE::Obj<RealSection>& slipRateSection =
- _fields->get("slip rate").section();
- assert(!slipRateSection.isNull());
+ double_array slipRateVertex(spaceDim);
+ const ALE::Obj<RealSection>& velRelSection =
+ _fields->get("relative velocity").section();
+ assert(!velRelSection.isNull());
const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
@@ -260,7 +414,7 @@
assert(!dispTIncrSection.isNull());
const ALE::Obj<RealSection>& orientationSection =
- fields->get("orientation").section();
+ _fields->get("orientation").section();
assert(!orientationSection.isNull());
const int numVertices = _cohesiveVertices.size();
@@ -270,13 +424,15 @@
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
- // Get slip
- assert(spaceDim == slipSection->getFiberDimension(v_fault));
- const double* slipVertex = slipSection->restrictPoint(v_fault);
+ // Get relative displacement
+ assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
+ const double* dispRelVertex = dispRelSection->restrictPoint(v_fault);
+ assert(dispRelVertex);
- // Get slip rate
- assert(spaceDim == slipRateSection->getFiberDimension(v_fault));
- const double* slipRateVertex = slipRateSection->restrictPoint(v_fault);
+ // Get relative velocity
+ assert(spaceDim == velRelSection->getFiberDimension(v_fault));
+ const double* velRelVertex = velRelSection->restrictPoint(v_fault);
+ assert(velRelVertex);
// Get orientation
assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
@@ -290,14 +446,19 @@
const double* lagrangeTIncrVertex =
dispTIncrSection->restrictPoint(v_lagrange);
- // Compute fault traction (Lagrange multiplier) at time t+dt in
- // fault coordinate system.
+ // Compute slip, slip rate, and fault traction (Lagrange
+ // multiplier) at time t+dt in fault coordinate system.
+ slipVertex = 0.0;
+ slipRateVertex = 0.0;
tractionTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
- tractionTpdtVertex[iDim] +=
- (lagrangeTVertex[jDim]+lagrangeTIncrVertex[jDim]) *
- orientationVertex[jDim*spaceDim+iDim];
+ slipVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
+ dispRelVertex[jDim];
+ slipRateVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
+ velRelVertex[jDim];
+ tractionTpdtVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
+ (lagrangeTVertex[jDim]+lagrangeTIncrVertex[jDim]);
} // for
} // for
@@ -359,7 +520,7 @@
assert(0 != _fields);
assert(0 != _friction);
- _updateSlipRate(*fields);
+ _updateVelRel(*fields);
_sensitivitySetup(jacobian);
// Update time step in friction (can vary).
@@ -371,22 +532,22 @@
double_array tractionTpdtVertex(spaceDim);
// Get sections
+ double_array dDispRelVertex(spaceDim);
+ const ALE::Obj<RealSection>& dispRelSection =
+ _fields->get("relative disp").section();
+ assert(!dispRelSection.isNull());
+
double_array dSlipVertex(spaceDim);
double_array slipVertex(spaceDim);
- const ALE::Obj<RealSection>& slipSection = _fields->get("slip").section();
- assert(!slipSection.isNull());
-
double_array slipRateVertex(spaceDim);
- const ALE::Obj<RealSection>& slipRateSection =
- _fields->get("slip rate").section();
- assert(!slipRateSection.isNull());
+ const ALE::Obj<RealSection>& velRelSection =
+ _fields->get("relative velocity").section();
+ assert(!velRelSection.isNull());
const ALE::Obj<RealSection>& orientationSection =
_fields->get("orientation").section();
assert(!orientationSection.isNull());
- double_array dispTVertexN(spaceDim);
- double_array dispTVertexP(spaceDim);
const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
@@ -423,24 +584,20 @@
} // switch
-#if 0 // DEBUGGING
- slipSection->view("SLIP");
- slipRateSection->view("SLIP RATE");
- //dispTSection->view("DISP (t)");
- //dispTIncrSection->view("DISP INCR (t->t+dt)");
-#endif
-
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
- // Get slip
- slipSection->restrictPoint(v_fault, &slipVertex[0], slipVertex.size());
+ // Get relative displacement
+ assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
+ const double* dispRelVertex = dispRelSection->restrictPoint(v_fault);
+ assert(dispRelVertex);
- // Get slip rate
- slipRateSection->restrictPoint(v_fault, &slipRateVertex[0],
- slipRateVertex.size());
+ // Get relative velocity
+ assert(spaceDim == velRelSection->getFiberDimension(v_fault));
+ const double* velRelVertex = velRelSection->restrictPoint(v_fault);
+ assert(velRelVertex);
// Get orientation
assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
@@ -455,10 +612,17 @@
const double* lagrangeTIncrVertex =
dispTIncrSection->restrictPoint(v_lagrange);
- // Compute Lagrange multiplier at time t+dt in fault coordinate system.
+ // Compute slip, slip rate, and Lagrange multiplier at time t+dt
+ // in fault coordinate system.
+ slipVertex = 0.0;
+ slipRateVertex = 0.0;
tractionTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
+ slipVertex[jDim] += orientationVertex[jDim*spaceDim+iDim] *
+ dispRelVertex[jDim];
+ slipRateVertex[jDim] += orientationVertex[jDim*spaceDim+iDim] *
+ velRelVertex[jDim];
tractionTpdtVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
(lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim]);
} // for
@@ -475,7 +639,7 @@
slipVertex, slipRateVertex,
tractionTpdtVertex);
- // Rotate traction to global coordinate system.
+ // Rotate traction back to global coordinate system.
dLagrangeTpdtVertexGlobal = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
@@ -505,30 +669,41 @@
// Update slip field based on solution of sensitivity problem and
// increment in Lagrange multipliers.
- const ALE::Obj<RealSection>& dispRelSection =
- _fields->get("sensitivity dispRel").section();
+ const ALE::Obj<RealSection>& sensDispRelSection =
+ _fields->get("sensitivity relative disp").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 slip
- slipSection->restrictPoint(v_fault, &slipVertex[0], slipVertex.size());
-
- // Get change in relative displacement.
+ // Get current relative displacement for updating.
+ assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
const double* dispRelVertex = dispRelSection->restrictPoint(v_fault);
assert(0 != dispRelVertex);
- assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
+ // Get change in relative displacement from sensitivity solve.
+ assert(spaceDim == sensDispRelSection->getFiberDimension(v_fault));
+ const double* sensDispRelVertex =
+ sensDispRelSection->restrictPoint(v_fault);
+ assert(sensDispRelVertex);
+
+ // Get orientation.
+ assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+ const double* orientationVertex =
+ orientationSection->restrictPoint(v_fault);
+ assert(orientationVertex);
+
// Get Lagrange multiplier at time t
assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
const double* lagrangeTVertex = dispTSection->restrictPoint(v_lagrange);
+ assert(lagrangeTVertex);
// Get Lagrange multiplier increment at time t
assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
const double* lagrangeTIncrVertex =
dispTIncrSection->restrictPoint(v_lagrange);
+ assert(lagrangeTIncrVertex);
// Get change in Lagrange multiplier.
dLagrangeTpdtSection->restrictPoint(v_fault, &dLagrangeTpdtVertex[0],
@@ -541,25 +716,47 @@
if (0.0 == dLagrangeMag)
continue; // No change, so continue
- // Compute change in slip.
- for (int iDim=0; iDim < spaceDim; ++iDim)
- dSlipVertex[iDim] = 2.0*dispRelVertex[iDim];
+ // Compute slip and change in slip in fault coordinates.
+ dSlipVertex = 0.0;
+ slipVertex = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
+ slipVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
+ dispRelVertex[jDim];
+ dSlipVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
+ 2.0*sensDispRelVertex[jDim];
+ } // for
+ } // for
+ // Compute normal traction in fault coordinates.
+ double tractionNormal = 0.0;
+ const int indexN = spaceDim - 1;
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
+ tractionNormal += orientationVertex[jDim*spaceDim+indexN] *
+ (lagrangeTVertex[indexN] + lagrangeTIncrVertex[indexN] +
+ dLagrangeTpdtVertex[indexN]);
+ } // for
+
// Do not allow fault interpenetration and set fault opening to
// zero if fault is under compression.
- // :TODO: FIX THIS
- const int indexN = spaceDim - 1;
- const double lagrangeTpdtNormal = lagrangeTVertex[indexN] +
- lagrangeTIncrVertex[indexN] + dLagrangeTpdtVertex[indexN];
- if (lagrangeTpdtNormal < -_zeroTolerance ||
+ if (tractionNormal < -_zeroTolerance ||
slipVertex[indexN] + dSlipVertex[indexN] < 0.0) {
dSlipVertex[indexN] = -slipVertex[indexN];
} // if
+ // Compute change relative displacement from change in slip.
+ dDispRelVertex = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
+ dDispRelVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ dSlipVertex[jDim];
+ } // for
+ } // for
+
// Set change in slip.
- assert(dSlipVertex.size() ==
- slipSection->getFiberDimension(v_fault));
- slipSection->updateAddPoint(v_fault, &dSlipVertex[0]);
+ assert(dDispRelVertex.size() ==
+ dispRelSection->getFiberDimension(v_fault));
+ dispRelSection->updateAddPoint(v_fault, &dDispRelVertex[0]);
// Update Lagrange multiplier increment.
assert(dLagrangeTpdtVertex.size() ==
@@ -567,7 +764,7 @@
dispTIncrSection->updateAddPoint(v_lagrange, &dLagrangeTpdtVertex[0]);
// Update displacement field
- dispTIncrVertexN = -0.5*dSlipVertex;
+ dispTIncrVertexN = -0.5*dDispRelVertex;
assert(dispTIncrVertexN.size() ==
dispTIncrSection->getFiberDimension(v_negative));
dispTIncrSection->updateAddPoint(v_negative, &dispTIncrVertexN[0]);
@@ -582,8 +779,8 @@
#if 0 // DEBUGGING
dLagrangeTpdtSection->view("AFTER dLagrange");
//dispTIncrSection->view("AFTER DISP INCR (t->t+dt)");
- slipSection->view("AFTER SLIP");
- slipRateSection->view("AFTER SLIP RATE");
+ dispRelSection->view("AFTER RELATIVE DISPLACEMENT");
+ velRelSection->view("AFTER RELATIVE VELOCITY");
#endif
} // constrainSolnSpace
@@ -595,6 +792,7 @@
topology::SolutionFields* const fields,
const topology::Field<topology::Mesh>& jacobian)
{ // adjustSolnLumped
+#if 0
/// Member prototype for _constrainSolnSpaceXD()
typedef void (pylith::faults::FaultCohesiveDyn::*constrainSolnSpace_fn_type)
(double_array*,
@@ -972,6 +1170,7 @@
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
#endif
+#endif
} // adjustSolnLumped
// ----------------------------------------------------------------------
@@ -992,13 +1191,26 @@
double scale = 0.0;
int fiberDim = 0;
if (0 == strcasecmp("slip", name)) {
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- return slip;
+ const topology::Field<topology::SubMesh>& dispRel =
+ _fields->get("relative disp");
+ _allocateBufferVectorField();
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ buffer.copy(dispRel);
+ buffer.label("slip");
+ _globalToFault(&buffer);
+ return buffer;
} else if (0 == strcasecmp("slip_rate", name)) {
- const topology::Field<topology::SubMesh>& slipRate =
- _fields->get("slip rate");
- return slipRate;
+ const topology::Field<topology::SubMesh>& velRel =
+ _fields->get("relative velocity");
+ _allocateBufferVectorField();
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ buffer.copy(velRel);
+ buffer.label("slip_rate");
+ _globalToFault(&buffer);
+ return buffer;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
const ALE::Obj<RealSection>& orientationSection = _fields->get(
@@ -1050,7 +1262,10 @@
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer =
_fields->get("buffer (vector)");
- _getInitialTractions(&buffer);
+ topology::Field<topology::SubMesh>& tractions =
+ _fields->get("initial tractions");
+ buffer.copy(tractions);
+ _globalToFault(&buffer);
return buffer;
} else if (0 == strcasecmp("traction", name)) {
@@ -1072,11 +1287,13 @@
throw std::runtime_error(msg.str());
} // else
+ // Should never get here.
+ throw std::logic_error("Unknown field in FaultCohesiveDyn::vertexField().");
+
// Satisfy return values
assert(0 != _fields);
const topology::Field<topology::SubMesh>& buffer = _fields->get(
"buffer (vector)");
- throw std::logic_error("Internal error.");
return buffer;
} // vertexField
@@ -1086,7 +1303,6 @@
pylith::faults::FaultCohesiveDyn::_setupInitialTractions(void)
{ // _setupInitialTractions
assert(0 != _normalizer);
- assert(0 != _quadrature);
// If no initial tractions specified, leave method
if (0 == _dbInitialTract)
@@ -1096,29 +1312,38 @@
const double pressureScale = _normalizer->pressureScale();
const double lengthScale = _normalizer->lengthScale();
- // Get quadrature information
- const int numQuadPts = _quadrature->numQuadPts();
- const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
- const double_array& quadWts = _quadrature->quadWts();
- assert(quadWts.size() == numQuadPts);
- double_array quadPtsGlobal(numQuadPts*spaceDim);
-
// Create section to hold initial tractions.
- _fields->add("initial forces", "initial_forces");
- topology::Field<topology::SubMesh>& forcesInitial =
- _fields->get("initial forces");
- topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- forcesInitial.cloneSection(slip);
- forcesInitial.scale(pressureScale);
- const ALE::Obj<RealSection>& forcesInitialSection = forcesInitial.section();
- assert(!forcesInitialSection.isNull());
- double_array forcesInitialCell(numBasis*spaceDim);
- double_array tractionQuadPt(spaceDim);
- UpdateAddVisitor forcesInitialVisitor(*forcesInitialSection,
- &forcesInitialCell[0]);
+ _fields->add("initial tractions", "initial_tractions");
+ topology::Field<topology::SubMesh>& initialTractions =
+ _fields->get("initial tractions");
+ topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
+ initialTractions.cloneSection(dispRel);
+ initialTractions.scale(pressureScale);
+ double_array initialTractionsVertex(spaceDim);
+ double_array initialTractionsVertexGlobal(spaceDim);
+ const ALE::Obj<RealSection>& initialTractionsSection =
+ initialTractions.section();
+ assert(!initialTractionsSection.isNull());
+
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
+ assert(!orientationSection.isNull());
+
+ const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
+ assert(0 != cs);
+
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+
+ double_array coordsVertex(spaceDim);
+ const ALE::Obj<RealSection>& coordsSection =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!coordsSection.isNull());
+
+
assert(0 != _dbInitialTract);
_dbInitialTract->open();
switch (spaceDim) { // switch
@@ -1143,124 +1368,61 @@
assert(0);
throw std::logic_error("Bad spatial dimension in Neumann.");
} // switch
-
- // Get cells associated with fault
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
- const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
- assert(0 != cs);
+ const int numVertices = _cohesiveVertices.size();
+ for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+ const int v_fault = _cohesiveVertices[iVertex].fault;
-#if !defined(PRECOMPUTE_GEOMETRY)
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
-#endif
+ coordsSection->restrictPoint(v_fault, &coordsVertex[0], coordsVertex.size());
- for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- // Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
-#else
- coordsVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
-#endif
+ assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+ const double* orientationVertex =
+ orientationSection->restrictPoint(v_fault);
+ assert(orientationVertex);
- const double_array& quadPtsNonDim = _quadrature->quadPts();
- quadPtsGlobal = quadPtsNonDim;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
- forcesInitialCell = 0.0;
+ _normalizer->dimensionalize(&coordsVertex[0], coordsVertex.size(),
+ lengthScale);
- // Loop over quadrature points in cell and query database
- for (int iQuadPt=0, index=0;
- iQuadPt < numQuadPts;
- ++iQuadPt, index+=spaceDim) {
+ initialTractionsVertex = 0.0;
+ int err = _dbInitialTract->query(&initialTractionsVertex[0],
+ initialTractionsVertex.size(),
+ &coordsVertex[0], coordsVertex.size(), cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find parameters for physical properties at \n" << "(";
+ for (int i = 0; i < spaceDim; ++i)
+ msg << " " << coordsVertex[i];
+ msg << ") in friction model " << label() << "\n"
+ << "using spatial database '" << _dbInitialTract->label() << "'.";
+ throw std::runtime_error(msg.str());
+ } // if
+ _normalizer->nondimensionalize(&initialTractionsVertex[0],
+ initialTractionsVertex.size(),
+ pressureScale);
- tractionQuadPt = 0.0;
- int err = _dbInitialTract->query(&tractionQuadPt[0], spaceDim,
- &quadPtsGlobal[index], spaceDim, cs);
- if (err) {
- std::ostringstream msg;
- msg << "Could not find parameters for physical properties at \n" << "(";
- for (int i = 0; i < spaceDim; ++i)
- msg << " " << quadPtsGlobal[index + i];
- msg << ") in friction model " << label() << "\n"
- << "using spatial database '" << _dbInitialTract->label() << "'.";
- throw std::runtime_error(msg.str());
- } // if
- tractionQuadPt /= pressureScale;
-
- // Get cell geometry information that depends on cell
- const double_array& basis = _quadrature->basis();
- const double_array& jacobianDet = _quadrature->jacobianDet();
-
- // Integrate tractions over cell.
- const double wt = quadWts[iQuadPt] * jacobianDet[iQuadPt];
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const double valI = wt*basis[iQuadPt*numBasis+iBasis];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const double valIJ = valI * basis[iQuadPt*numBasis+jBasis];
- for (int iDim=0; iDim < spaceDim; ++iDim)
- forcesInitialCell[iBasis*spaceDim+iDim] +=
- tractionQuadPt[iDim] * valIJ;
- } // for
+ // Rotate tractions from fault coordinate system to global
+ // coordinate system
+ initialTractionsVertexGlobal = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
+ initialTractionsVertexGlobal[iDim] +=
+ orientationVertex[jDim*spaceDim+iDim] *
+ initialTractionsVertex[jDim];
} // for
} // for
- // Assemble cell contribution into field
- forcesInitialVisitor.clear();
- faultSieveMesh->updateClosure(*c_iter, forcesInitialVisitor);
+
+ assert(initialTractionsVertexGlobal.size() ==
+ initialTractionsSection->getFiberDimension(v_fault));
+ initialTractionsSection->updatePoint(v_fault,
+ &initialTractionsVertexGlobal[0]);
} // for
+
// Close properties database
_dbInitialTract->close();
- forcesInitial.complete(); // Assemble contributions
+ initialTractions.complete(); // Assemble contributions
- // Rotate forces from fault coordinate system to global coordinate system
- const int orientationSize = spaceDim * spaceDim;
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
-
- double_array forcesInitialVertexFault(spaceDim);
- double_array forcesInitialVertexGlobal(spaceDim);
-
- const int numVertices = _cohesiveVertices.size();
- for (int iVertex=0; iVertex < numVertices; ++iVertex) {
- const int v_fault = _cohesiveVertices[iVertex].fault;
-
- assert(orientationSize == orientationSection->getFiberDimension(v_fault));
- assert(spaceDim == forcesInitialSection->getFiberDimension(v_fault));
-
- const double* orientationVertex =
- orientationSection->restrictPoint(v_fault);
-
- forcesInitialSection->restrictPoint(v_fault,
- &forcesInitialVertexFault[0],
- forcesInitialVertexFault.size());
-
- forcesInitialVertexGlobal = 0.0;
- for (int iDim = 0; iDim < spaceDim; ++iDim)
- for (int kDim = 0; kDim < spaceDim; ++kDim)
- forcesInitialVertexGlobal[iDim] +=
- forcesInitialVertexFault[kDim] *
- orientationVertex[kDim*spaceDim+iDim];
-
- assert(forcesInitialVertexGlobal.size() ==
- forcesInitialSection->getFiberDimension(v_fault));
- forcesInitialSection->updatePoint(v_fault, &forcesInitialVertexGlobal[0]);
- } // for
-
- //forcesInitial.view("INITIAL FORCES"); // DEBUGGING
+ //initalTractions.view("INITIAL FORCES"); // DEBUGGING
} // _setupInitialTractions
// ----------------------------------------------------------------------
@@ -1276,82 +1438,16 @@
assert(0 != _normalizer);
// Fiber dimension of tractions matches spatial dimension.
- const int fiberDim = _quadrature->spaceDim();
- double_array tractionsVertex(fiberDim);
+ const int spaceDim = _quadrature->spaceDim();
+ double_array tractionsVertex(spaceDim);
// Get sections.
const ALE::Obj<RealSection>& dispTSection = dispT.section();
assert(!dispTSection.isNull());
- // Allocate buffer for tractions field (if necessary).
- const ALE::Obj<RealSection>& tractionsSection = tractions->section();
- if (tractionsSection.isNull()) {
- ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- //logger.stagePush("Fault");
-
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- tractions->cloneSection(slip);
-
- //logger.stagePop();
- } // if
- const double pressureScale = _normalizer->pressureScale();
- tractions->label("traction");
- tractions->scale(pressureScale);
- tractions->zero();
-
- const int numVertices = _cohesiveVertices.size();
- for (int iVertex=0; iVertex < numVertices; ++iVertex) {
- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
- const int v_fault = _cohesiveVertices[iVertex].fault;
-
- assert(fiberDim == dispTSection->getFiberDimension(v_lagrange));
- assert(fiberDim == tractionsSection->getFiberDimension(v_fault));
-
- const double* dispTVertex = dispTSection->restrictPoint(v_lagrange);
- assert(0 != dispTVertex);
-
- // :TODO: Rotate to fault coordinate system
-#if 0 // :TODO: FIX THIS
- for (int i=0; i < fiberDim; ++i)
- tractionsVertex[i] = dispTVertex[i] / areaVertex[0];
-#endif
-
- assert(tractionsVertex.size() ==
- tractionsSection->getFiberDimension(v_fault));
- tractionsSection->updatePoint(v_fault, &tractionsVertex[0]);
- } // for
-
- PetscLogFlops(numVertices * (1 + fiberDim) );
-
-#if 0 // DEBUGGING
- tractions->view("TRACTIONS");
-#endif
-
-} // _calcTractions
-
-// ----------------------------------------------------------------------
-// Compute initial tractions on fault surface.
-void
-pylith::faults::FaultCohesiveDyn::_getInitialTractions(
- topology::Field<topology::SubMesh>* tractions)
-{ // _getInitialTractions
- assert(0 != tractions);
- assert(0 != _faultMesh);
- assert(0 != _fields);
- assert(0 != _normalizer);
-
- // Fiber dimension of tractions matches spatial dimension.
- const int spaceDim = _quadrature->spaceDim();
- double_array tractionsVertexGlobal(spaceDim);
- double_array tractionsVertexFault(spaceDim);
-
- // Get sections.
const ALE::Obj<RealSection>& orientationSection =
_fields->get("orientation").section();
assert(!orientationSection.isNull());
- const ALE::Obj<RealSection>& forcesInitialSection =
- _fields->get("initial forces").section();
- assert(!forcesInitialSection.isNull());
// Allocate buffer for tractions field (if necessary).
const ALE::Obj<RealSection>& tractionsSection = tractions->section();
@@ -1359,81 +1455,73 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
//logger.stagePush("Fault");
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- tractions->cloneSection(slip);
+ const topology::Field<topology::SubMesh>& dispRel =
+ _fields->get("relative disp");
+ tractions->cloneSection(dispRel);
//logger.stagePop();
} // if
const double pressureScale = _normalizer->pressureScale();
- tractions->label("initial_traction");
+ tractions->label("traction");
tractions->scale(pressureScale);
tractions->zero();
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+ const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
- assert(spaceDim == forcesInitialSection->getFiberDimension(v_fault));
- assert(spaceDim == tractionsSection->getFiberDimension(v_fault));
+ assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
+ const double* dispTVertex = dispTSection->restrictPoint(v_lagrange);
+ assert(dispTVertex);
- const double* forcesInitialVertex =
- forcesInitialSection->restrictPoint(v_fault);
- assert(0 != forcesInitialVertex);
+ assert(spaceDim*spaceDim ==
+ orientationSection->getFiberDimension(v_fault));
const double* orientationVertex =
orientationSection->restrictPoint(v_fault);
- assert(0 != orientationVertex);
+ assert(orientationVertex);
-#if 0 // :TODO: FIX THIS
- for (int i = 0; i < spaceDim; ++i)
- tractionsVertexGlobal[i] = forcesInitialVertex[i] / areaVertex[0];
-#endif
+ // Rotate tractions to fault coordinate system.
+ tractionsVertex = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
+ tractionsVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
+ dispTVertex[jDim];
+ } // for
+ } // for
- // Rotate from global coordinate system to fault coordinate system
- tractionsVertexFault = 0.0;
- for (int iDim = 0; iDim < spaceDim; ++iDim)
- for (int kDim = 0; kDim < spaceDim; ++kDim)
- tractionsVertexFault[iDim] +=
- tractionsVertexGlobal[kDim] * orientationVertex[iDim*spaceDim+kDim];
-
- assert(tractionsVertexFault.size() ==
+ assert(tractionsVertex.size() ==
tractionsSection->getFiberDimension(v_fault));
- tractionsSection->updatePoint(v_fault, &tractionsVertexFault[0]);
+ tractionsSection->updatePoint(v_fault, &tractionsVertex[0]);
} // for
PetscLogFlops(numVertices * (1 + spaceDim) );
#if 0 // DEBUGGING
- tractions->view("INITIAL TRACTIONS");
+ tractions->view("TRACTIONS");
#endif
-} // _getInitialTractions
+} // _calcTractions
// ----------------------------------------------------------------------
// Update slip rate associated with Lagrange vertex k corresponding
// to diffential velocity between conventional vertices i and j.
void
-pylith::faults::FaultCohesiveDyn::_updateSlipRate(const topology::SolutionFields& fields)
-{ // _updateSlipRate
+pylith::faults::FaultCohesiveDyn::_updateVelRel(const topology::SolutionFields& fields)
+{ // _updateVelRel
assert(0 != _fields);
const int spaceDim = _quadrature->spaceDim();
// Get section information
- double_array velocityVertexN(spaceDim);
- double_array velocityVertexP(spaceDim);
const ALE::Obj<RealSection>& velocitySection =
fields.get("velocity(t)").section();
assert(!velocitySection.isNull());
- double_array slipRateVertex(spaceDim);
- const ALE::Obj<RealSection>& slipRateSection =
- _fields->get("slip rate").section();
+ double_array velRelVertex(spaceDim);
+ const ALE::Obj<RealSection>& velRelSection =
+ _fields->get("relative velocity").section();
- double_array orientationVertex(spaceDim*spaceDim);
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
- assert(!orientationSection.isNull());
-
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
@@ -1450,36 +1538,20 @@
assert(0 != velocityVertexP);
assert(spaceDim == velocitySection->getFiberDimension(v_positive));
- const double* orientationVertex = orientationSection->restrictPoint(v_fault);
- assert(0 != orientationVertex);
- assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+ // Compute relative velocity
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const double value = velocityVertexP[iDim] - velocityVertexN[iDim];
+ velRelVertex[iDim] = fabs(value) > _zeroTolerance ? value : 0.0;
+ } // for
- slipRateVertex = 0.0;
- // Velocity for negative vertex.
- for (int iDim = 0; iDim < spaceDim; ++iDim)
- for (int kDim = 0; kDim < spaceDim; ++kDim)
- slipRateVertex[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[iDim*spaceDim+kDim];
-
- // Limit velocity to resolvable range
- for (int iDim = 0; iDim < spaceDim; ++iDim)
- if (fabs(slipRateVertex[iDim]) < _zeroTolerance)
- slipRateVertex[iDim] = 0.0;
-
// Update slip rate field.
- assert(slipRateVertex.size() ==
- slipRateSection->getFiberDimension(v_fault));
- slipRateSection->updatePoint(v_fault, &slipRateVertex[0]);
+ assert(velRelVertex.size() ==
+ velRelSection->getFiberDimension(v_fault));
+ velRelSection->updatePoint(v_fault, &velRelVertex[0]);
} // for
PetscLogFlops(numVertices*spaceDim*spaceDim*4);
-} // _updateSlipRate
+} // _updateVelRel
// ----------------------------------------------------------------------
// Setup sensitivity problem to compute change in slip given change in Lagrange multipliers.
@@ -1496,9 +1568,9 @@
_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);
+ const topology::Field<topology::SubMesh>& dispRel =
+ _fields->get("relative disp");
+ solution.cloneSection(dispRel);
solution.createScatter(solution.mesh());
} // if
const topology::Field<topology::SubMesh>& solution =
@@ -1513,13 +1585,13 @@
} // if
if (!_fields->hasField("sensitivity dispRel")) {
- _fields->add("sensitivity dispRel", "sensitivity_disprel");
+ _fields->add("sensitivity relative disp", "sensitivity_relative_disp");
topology::Field<topology::SubMesh>& dispRel =
- _fields->get("sensitivity dispRel");
+ _fields->get("sensitivity relative disp");
dispRel.cloneSection(solution);
} // if
topology::Field<topology::SubMesh>& dispRel =
- _fields->get("sensitivity dispRel");
+ _fields->get("sensitivity relative disp");
dispRel.zero();
if (!_fields->hasField("sensitivity dLagrange")) {
@@ -1783,7 +1855,7 @@
const ALE::Obj<RealSection>& solutionSection =
_fields->get("sensitivity solution").section();
const ALE::Obj<RealSection>& dispRelSection =
- _fields->get("sensitivity dispRel").section();
+ _fields->get("sensitivity relative disp").section();
const double sign = (negativeSide) ? -1.0 : 1.0;
Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.hh 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.hh 2011-09-30 23:02:10 UTC (rev 18998)
@@ -146,18 +146,13 @@
void _calcTractions(topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& solution);
- /** Get initial tractions on fault surface.
+ /** Update relative velocity associated with Lagrange vertex k
+ * corresponding to diffential velocity between conventional
+ * vertices i and j.
*
- * @param tractions Field for tractions.
- */
- void _getInitialTractions(topology::Field<topology::SubMesh>* tractions);
-
- /** Update slip rate associated with Lagrange vertex k corresponding
- * to diffential velocity between conventional vertices i and j.
- *
* @param fields Solution fields.
*/
- void _updateSlipRate(const topology::SolutionFields& fields);
+ void _updateVelRel(const topology::SolutionFields& fields);
/** Setup sensitivity problem to compute change in slip given change
* in Lagrange multipliers.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveKin.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveKin.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -129,20 +129,20 @@
const int setupEvent = _logger->eventId("FaIR setup");
_logger->eventBegin(setupEvent);
- topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- slip.zero();
+ topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
+ dispRel.zero();
// Compute slip field at current time step
const srcs_type::const_iterator srcsEnd = _eqSrcs.end();
for (srcs_type::iterator s_iter = _eqSrcs.begin(); s_iter != srcsEnd; ++s_iter) {
EqKinSrc* src = s_iter->second;
assert(0 != src);
if (t >= src->originTime())
- src->slip(&slip, t);
+ src->slip(&dispRel, t);
} // for
- // Transform slip from local (fault) coordinate system to global
- // coordinate system
- _slipFaultToGlobal();
+ // Transform slip from local (fault) coordinate system to relative
+ // displacement field in global coordinate system
+ _faultToGlobal(&dispRel);
_logger->eventEnd(setupEvent);
@@ -170,8 +170,15 @@
double scale = 0.0;
int fiberDim = 0;
if (0 == strcasecmp("slip", name)) {
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- return slip;
+ const topology::Field<topology::SubMesh>& dispRel =
+ _fields->get("relative disp");
+ _allocateBufferVectorField();
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ buffer.copy(dispRel);
+ buffer.label("slip");
+ _globalToFault(&buffer);
+ return buffer;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
const ALE::Obj<RealSection>& orientationSection = _fields->get(
@@ -272,6 +279,9 @@
} // else
+ // Should never get here.
+ throw std::logic_error("Unknown field in FaultCohesiveKin::vertexField().");
+
// Satisfy return values
assert(0 != _fields);
const topology::Field<topology::SubMesh>& buffer = _fields->get(
Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -107,18 +107,18 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
//logger.stagePush("Fault");
- // Allocate slip field
+ // Allocate dispRel field
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
faultSieveMesh->depthStratum(0);
assert(!vertices.isNull());
- _fields->add("slip", "slip");
- topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- slip.newSection(vertices, cs->spaceDim());
- slip.allocate();
- slip.vectorFieldType(topology::FieldBase::VECTOR);
- slip.scale(_normalizer->lengthScale());
+ _fields->add("relative disp", "relative_disp");
+ topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
+ dispRel.newSection(vertices, cs->spaceDim());
+ dispRel.allocate();
+ dispRel.vectorFieldType(topology::FieldBase::VECTOR);
+ dispRel.scale(_normalizer->lengthScale());
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
faultSieveMesh->heightStratum(0);
@@ -267,11 +267,12 @@
RestrictVisitor coordsVisitor(*coordinates,
coordinatesCell.size(), &coordinatesCell[0]);
- double_array slipCell(numBasis*spaceDim);
- topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- const ALE::Obj<RealSection>& slipSection = slip.section();
- assert(!slipSection.isNull());
- RestrictVisitor slipVisitor(*slipSection, slipCell.size(), &slipCell[0]);
+ double_array dispRelCell(numBasis*spaceDim);
+ topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
+ const ALE::Obj<RealSection>& dispRelSection = dispRel.section();
+ assert(!dispRelSection.isNull());
+ RestrictVisitor dispRelVisitor(*dispRelSection,
+ dispRelCell.size(), &dispRelCell[0]);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -308,8 +309,8 @@
dispTIncrVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
- slipVisitor.clear();
- faultSieveMesh->restrictClosure(c_fault, slipVisitor);
+ dispRelVisitor.clear();
+ faultSieveMesh->restrictClosure(c_fault, dispRelVisitor);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -370,7 +371,7 @@
// Lagrange constraint
residualCell[iBL + iDim] += valIJ *
(dispTpdtCell[jBP + iDim] - dispTpdtCell[jBN + iDim] -
- slipCell[jBasis*spaceDim+iDim]);
+ dispRelCell[jBasis*spaceDim+iDim]);
#if 0
std::cout << "iBasis: " << iBasis
@@ -378,7 +379,7 @@
<< ", iDim: " << iDim
<< ", valIJ: " << valIJ
<< ", jacobianDet: " << jacobianDet[iQuad]
- << ", slip: " << slipCell[jBasis*spaceDim+iDim]
+ << ", dispRel: " << dispRelCell[jBasis*spaceDim+iDim]
<< ", dispP: " << dispTpdtCell[jBP + iDim]
<< ", dispN: " << dispTpdtCell[jBN + iDim]
<< ", dispL: " << dispTpdtCell[jBL + iDim]
@@ -1515,8 +1516,9 @@
// Allocate orientation field.
_fields->add("orientation", "orientation");
topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- orientation.newSection(slip, orientationSize);
+ const topology::Field<topology::SubMesh>& dispRel =
+ _fields->get("relative disp");
+ orientation.newSection(dispRel, orientationSize);
const ALE::Obj<RealSection>& orientationSection = orientation.section();
assert(!orientationSection.isNull());
// Create subspaces for along-strike, up-dip, and normal directions
@@ -1760,8 +1762,9 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
//logger.stagePush("Fault");
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- tractions->cloneSection(slip);
+ const topology::Field<topology::SubMesh>& dispRel =
+ _fields->get("relative disp");
+ tractions->cloneSection(dispRel);
//logger.stagePop();
} // if
@@ -1780,13 +1783,15 @@
assert(0 != dispTVertex);
assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
- const double* orientationVertex = orientationSection->restrictPoint(v_fault);
+ const double* orientationVertex =
+ orientationSection->restrictPoint(v_fault);
assert(orientationVertex);
tractionsVertex = 0.0;
- for (int i=0; i < spaceDim; ++i)
- for (int j=0; j < spaceDim; ++j)
- tractionsVertex[i] += orientationVertex[i*spaceDim+j] * dispTVertex[j];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < spaceDim; ++jDim)
+ tractionsVertex[iDim] +=
+ orientationVertex[iDim*spaceDim+jDim] * dispTVertex[jDim];
assert(tractionsVertex.size() ==
tractionsSection->getFiberDimension(v_fault));
@@ -1801,22 +1806,22 @@
} // _calcTractionsChange
// ----------------------------------------------------------------------
-// Transform slip field from local (fault) coordinate system to
+// Transform field from local (fault) coordinate system to
// global coordinate system.
void
-pylith::faults::FaultCohesiveLagrange::_slipFaultToGlobal(void)
-{ // _slipFaultToGlobal
+pylith::faults::FaultCohesiveLagrange::_faultToGlobal(topology::Field<topology::SubMesh>* field)
+{ // _faultToGlobal
+ assert(field);
assert(0 != _faultMesh);
assert(0 != _fields);
- // Fiber dimension of tractions matches spatial dimension.
+ // Fiber dimension of vector field matches spatial dimension.
const int spaceDim = _quadrature->spaceDim();
- double_array slipGlobalVertex(spaceDim);
+ double_array fieldVertexGlobal(spaceDim);
// Get sections.
- const ALE::Obj<RealSection>& slipSection =
- _fields->get("slip").section();
- assert(!slipSection.isNull());
+ const ALE::Obj<RealSection>& fieldSection = field->section();
+ assert(!fieldSection.isNull());
const ALE::Obj<RealSection>& orientationSection =
_fields->get("orientation").section();
@@ -1826,33 +1831,85 @@
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_fault = _cohesiveVertices[iVertex].fault;
- assert(spaceDim == slipSection->getFiberDimension(v_fault));
- const double* slipFaultVertex = slipSection->restrictPoint(v_fault);
- assert(slipFaultVertex);
+ assert(spaceDim == fieldSection->getFiberDimension(v_fault));
+ const double* fieldVertexFault = fieldSection->restrictPoint(v_fault);
+ assert(fieldVertexFault);
assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
const double* orientationVertex = orientationSection->restrictPoint(v_fault);
assert(orientationVertex);
- slipGlobalVertex = 0.0;
- for (int i=0; i < spaceDim; ++i)
- for (int j=0; j < spaceDim; ++j)
- slipGlobalVertex[i] +=
- orientationVertex[j*spaceDim+i] * slipFaultVertex[j];
+ fieldVertexGlobal = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < spaceDim; ++jDim)
+ fieldVertexGlobal[iDim] +=
+ orientationVertex[jDim*spaceDim+iDim] * fieldVertexFault[jDim];
- assert(slipGlobalVertex.size() ==
- slipSection->getFiberDimension(v_fault));
- slipSection->updatePoint(v_fault, &slipGlobalVertex[0]);
+ assert(fieldVertexGlobal.size() ==
+ fieldSection->getFiberDimension(v_fault));
+ fieldSection->updatePoint(v_fault, &fieldVertexGlobal[0]);
} // for
PetscLogFlops(numVertices * (2*spaceDim*spaceDim) );
#if 0 // DEBUGGING
- slipSection->view("SLIP (GLOBAL)");
+ field->view("FIELD (GLOBAL)");
#endif
-} // _slipFaultToGlobal
+} // _faultToGlobal
// ----------------------------------------------------------------------
+// Transform field from global coordinate system to local (fault)
+// coordinate system.
+void
+pylith::faults::FaultCohesiveLagrange::_globalToFault(topology::Field<topology::SubMesh>* field)
+{ // _globalToFault
+ assert(field);
+ assert(0 != _faultMesh);
+ assert(0 != _fields);
+
+ // Fiber dimension of vector field matches spatial dimension.
+ const int spaceDim = _quadrature->spaceDim();
+ double_array fieldVertexFault(spaceDim);
+
+ // Get sections.
+ const ALE::Obj<RealSection>& fieldSection = field->section();
+ assert(!fieldSection.isNull());
+
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
+ assert(!orientationSection.isNull());
+
+ const int numVertices = _cohesiveVertices.size();
+ for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+ const int v_fault = _cohesiveVertices[iVertex].fault;
+
+ assert(spaceDim == fieldSection->getFiberDimension(v_fault));
+ const double* fieldVertexGlobal = fieldSection->restrictPoint(v_fault);
+ assert(fieldVertexGlobal);
+
+ assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+ const double* orientationVertex = orientationSection->restrictPoint(v_fault);
+ assert(orientationVertex);
+
+ fieldVertexFault = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < spaceDim; ++jDim)
+ fieldVertexFault[iDim] +=
+ orientationVertex[iDim*spaceDim+jDim] * fieldVertexGlobal[jDim];
+
+ assert(fieldVertexFault.size() ==
+ fieldSection->getFiberDimension(v_fault));
+ fieldSection->updatePoint(v_fault, &fieldVertexFault[0]);
+ } // for
+
+ PetscLogFlops(numVertices * (2*spaceDim*spaceDim) );
+
+#if 0 // DEBUGGING
+ field->view("FIELD (FAULT)");
+#endif
+} // _faultToGlobal
+
+// ----------------------------------------------------------------------
// Allocate buffer for vector field.
void
pylith::faults::FaultCohesiveLagrange::_allocateBufferVectorField(void)
@@ -1864,12 +1921,14 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Output");
- // Create vector field; use same shape/chart as slip field.
+ // Create vector field; use same shape/chart as relative
+ // displacement field.
assert(0 != _faultMesh);
_fields->add("buffer (vector)", "buffer");
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- buffer.cloneSection(slip);
+ const topology::Field<topology::SubMesh>& dispRel =
+ _fields->get("relative disp");
+ buffer.cloneSection(dispRel);
buffer.zero();
assert(buffer.vectorFieldType() == topology::FieldBase::VECTOR);
Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.hh 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.hh 2011-09-30 23:02:10 UTC (rev 18998)
@@ -234,11 +234,20 @@
void _calcTractionsChange(topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& solution);
- /** Transform slip field from local (fault) coordinate system to
+ /** Transform field from local (fault) coordinate system to
* global coordinate system.
+ *
+ * @param field Field to transform.
*/
- void _slipFaultToGlobal(void);
+ void _faultToGlobal(topology::Field<topology::SubMesh>* field);
+ /** Transform field from global coordinate system to local (fault)
+ * coordinate system.
+ *
+ * @param field Field to transform.
+ */
+ void _globalToFault(topology::Field<topology::SubMesh>* field);
+
/// Allocate buffer for vector field.
void _allocateBufferVectorField(void);
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -155,28 +155,28 @@
} // for
} // for
- // Initial forces/tractions
+ // Initial tractions
if (0 != fault._dbInitialTract) {
- //fault._fields->get("initial forces").view("INITIAL FORCES"); // DEBUGGING
- const ALE::Obj<RealSection>& forcesInitialSection =
- fault._fields->get("initial forces").section();
- CPPUNIT_ASSERT(!forcesInitialSection.isNull());
+ //fault._fields->get("initial tractions").view("INITIAL TRACTIONS"); // DEBUGGING
+ const ALE::Obj<RealSection>& initialTractionsSection =
+ fault._fields->get("initial tractions").section();
+ CPPUNIT_ASSERT(!initialTractionsSection.isNull());
const int spaceDim = _data->spaceDim;
iVertex = 0;
for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex) {
- const int fiberDim = forcesInitialSection->getFiberDimension(*v_iter);
+ const int fiberDim = initialTractionsSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const double* forcesInitialVertex =
- forcesInitialSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != forcesInitialVertex);
+ const double* initialTractionsVertex =
+ initialTractionsSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(initialTractionsVertex);
const double tolerance = 1.0e-06;
for (int i = 0; i < spaceDim; ++i) {
const int index = iVertex * spaceDim + i;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->forcesInitial[index],
- forcesInitialVertex[i], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->initialTractions[index],
+ initialTractionsVertex[i], tolerance);
} // for
} // for
} // if
@@ -252,8 +252,6 @@
// Get fault vertex info
const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- SieveSubMesh::renumbering_type& renumbering =
- faultSieveMesh->getRenumbering();
const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
faultSieveMesh->depthStratum(0);
CPPUNIT_ASSERT(!vertices.isNull());
@@ -366,8 +364,6 @@
// Get fault vertex info
const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- SieveSubMesh::renumbering_type& renumbering =
- faultSieveMesh->getRenumbering();
const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
faultSieveMesh->depthStratum(0);
CPPUNIT_ASSERT(!vertices.isNull());
@@ -490,8 +486,6 @@
// Get fault vertex info
const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- SieveSubMesh::renumbering_type& renumbering =
- faultSieveMesh->getRenumbering();
const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
faultSieveMesh->depthStratum(0);
CPPUNIT_ASSERT(!vertices.isNull());
@@ -629,22 +623,27 @@
int fiberDim = tractionsSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
const double* tractionsVertex = tractionsSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != tractionsVertex);
+ CPPUNIT_ASSERT(tractionsVertex);
- fiberDim = dispSection->getFiberDimension(meshVertex);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const double* dispVertex = dispSection->restrictPoint(meshVertex);
- CPPUNIT_ASSERT(0 != dispVertex);
+ const double* tractionsVertexGlobalE =
+ dispSection->restrictPoint(meshVertex);
+ CPPUNIT_ASSERT(tractionsVertexGlobalE);
+ const double* orientationVertex =
+ &_data->orientation[iVertex*spaceDim*spaceDim];
+ CPPUNIT_ASSERT(orientationVertex);
- const double scale = 1.0 / _data->area[iVertex];
for (int iDim=0; iDim < spaceDim; ++iDim) {
- const double tractionE = dispVertex[iDim] * scale;
+ double tractionE = 0.0;
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
+ tractionE +=
+ orientationVertex[jDim*spaceDim+iDim]*tractionsVertexGlobalE[jDim];
+ } // for
if (tractionE != 0.0)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, tractionsVertex[iDim]/tractionE,
tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionE, tractionsVertex[iDim],
- tolerance);
+ tolerance);
} // for
} // for
} // testCalcTractions
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynData.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynData.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynData.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -40,7 +40,7 @@
fieldIncrOpen(0),
jacobian(0),
orientation(0),
- forcesInitial(0),
+ initialTractions(0),
area(0),
fieldIncrSlipE(0),
slipSlipE(0),
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynData.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynData.hh 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynData.hh 2011-09-30 23:02:10 UTC (rev 18998)
@@ -75,7 +75,7 @@
//@{
double* orientation; ///< Expected values for fault orientation.
double* area; ///< Expected values for fault area.
- double* forcesInitial; ///< Expected values for initial forces.
+ double* initialTractions; ///< Expected values for initial tractions.
double* fieldIncrSlipE; ///< Expected values for solution increment for slipping case.
double* slipSlipE; ///< Expected values for slip for slipping case.
double* fieldIncrOpenE; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -1345,11 +1345,11 @@
1.0, 1.0, 1.0, 1.0
};
-const double pylith::faults::CohesiveDynDataHex8::_forcesInitial[] = {
- 3.063397471, -1.063397471, +2.063397471,
- 3.121132498, -1.121132498, +2.121132498,
- 3.178867525, -1.178867525, +2.178867525,
- 3.236602552, -1.236602552, +2.236602552,
+const double pylith::faults::CohesiveDynDataHex8::_initialTractions[] = {
+ 3.0, -1.0, +2.0,
+ 3.1, -1.1, +2.1,
+ 3.2, -1.2, +2.2,
+ 3.3, -1.3, +2.3,
};
@@ -1524,7 +1524,7 @@
jacobian = const_cast<double*>(_jacobian);
orientation = const_cast<double*>(_orientation);
area = const_cast<double*>(_area);
- forcesInitial = const_cast<double*>(_forcesInitial);
+ initialTractions = const_cast<double*>(_initialTractions);
constraintVertices = const_cast<int*>(_constraintVertices);
numConstraintVert = _numConstraintVert;
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.hh 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.hh 2011-09-30 23:02:10 UTC (rev 18998)
@@ -67,7 +67,7 @@
static const double _orientation[]; ///< Expected values for fault orientation.
static const double _area[]; ///< Expected values for fault area.
- static const double _forcesInitial[]; ///< Expected values for initial forces.
+ static const double _initialTractions[]; ///< Expected values for initial tractions.
static const double _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const double _slipSlipE[]; ///< Expected values for slip for slip case.
static const double _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -322,9 +322,9 @@
1.0,
};
-const double pylith::faults::CohesiveDynDataQuad4::_forcesInitial[] = {
- 2.05, -1.05,
- 2.05, -1.05,
+const double pylith::faults::CohesiveDynDataQuad4::_initialTractions[] = {
+ 2.0, -1.0,
+ 2.1, -1.1,
};
@@ -448,7 +448,7 @@
jacobian = const_cast<double*>(_jacobian);
orientation = const_cast<double*>(_orientation);
area = const_cast<double*>(_area);
- forcesInitial = const_cast<double*>(_forcesInitial);
+ initialTractions = const_cast<double*>(_initialTractions);
constraintVertices = const_cast<int*>(_constraintVertices);
numConstraintVert = _numConstraintVert;
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.hh 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.hh 2011-09-30 23:02:10 UTC (rev 18998)
@@ -67,7 +67,7 @@
static const double _orientation[]; ///< Expected values for fault orientation.
static const double _area[]; ///< Expected values for fault area.
- static const double _forcesInitial[]; ///< Expected values for initial forces.
+ static const double _initialTractions[]; ///< Expected values for initial tractions.
static const double _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const double _slipSlipE[]; ///< Expected values for slip for slip case.
static const double _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -476,10 +476,10 @@
1.0/3.0,
};
-const double pylith::faults::CohesiveDynDataTet4::_forcesInitial[] = {
- 3.1/3.0, -1.1/3.0, +2.1/3.0,
- 3.1/3.0, -1.1/3.0, +2.1/3.0,
- 3.1/3.0, -1.1/3.0, +2.1/3.0,
+const double pylith::faults::CohesiveDynDataTet4::_initialTractions[] = {
+ 3.0, -1.0, +2.0,
+ 3.1, -1.1, +2.1,
+ 3.2, -1.2, +2.2,
};
@@ -608,7 +608,7 @@
jacobian = const_cast<double*>(_jacobian);
orientation = const_cast<double*>(_orientation);
area = const_cast<double*>(_area);
- forcesInitial = const_cast<double*>(_forcesInitial);
+ initialTractions = const_cast<double*>(_initialTractions);
constraintVertices = const_cast<int*>(_constraintVertices);
numConstraintVert = _numConstraintVert;
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.hh 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.hh 2011-09-30 23:02:10 UTC (rev 18998)
@@ -67,7 +67,7 @@
static const double _orientation[]; ///< Expected values for fault orientation.
static const double _area[]; ///< Expected values for fault area.
- static const double _forcesInitial[]; ///< Expected values for initial forces.
+ static const double _initialTractions[]; ///< Expected values for initial tractions.
static const double _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const double _slipSlipE[]; ///< Expected values for slip for slip case.
static const double _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -252,9 +252,9 @@
1.0,
};
-const double pylith::faults::CohesiveDynDataTri3::_forcesInitial[] = {
- 2.05, -1.05,
- 2.05, -1.05,
+const double pylith::faults::CohesiveDynDataTri3::_initialTractions[] = {
+ 2.0, -1.0,
+ 2.1, -1.1,
};
@@ -366,7 +366,7 @@
jacobian = const_cast<double*>(_jacobian);
orientation = const_cast<double*>(_orientation);
area = const_cast<double*>(_area);
- forcesInitial = const_cast<double*>(_forcesInitial);
+ initialTractions = const_cast<double*>(_initialTractions);
constraintVertices = const_cast<int*>(_constraintVertices);
numConstraintVert = _numConstraintVert;
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.hh 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.hh 2011-09-30 23:02:10 UTC (rev 18998)
@@ -67,7 +67,7 @@
static const double _orientation[]; ///< Expected values for fault orientation.
static const double _area[]; ///< Expected values for fault area.
- static const double _forcesInitial[]; ///< Expected values for initial forces.
+ static const double _initialTractions[]; ///< Expected values for initial tractions.
static const double _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const double _slipSlipE[]; ///< Expected values for slip for slip case.
static const double _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc 2011-09-30 23:02:10 UTC (rev 18998)
@@ -438,10 +438,10 @@
13, 14, 15
};
-const double pylith::faults::CohesiveDynDataTri3d::_forcesInitial[] = {
- 3.15*1.4142135623730951, 1.00*1.41421356237309,
- 2.05, -1.05,
- 1.10, 2.10,
+const double pylith::faults::CohesiveDynDataTri3d::_initialTractions[] = {
+ 3.0*0.70710678118654757, 0.70710678118654757,
+ 2.1, -1.1,
+ 1.2, 2.2,
};
@@ -572,7 +572,7 @@
jacobian = const_cast<double*>(_jacobian);
orientation = const_cast<double*>(_orientation);
area = const_cast<double*>(_area);
- forcesInitial = const_cast<double*>(_forcesInitial);
+ initialTractions = const_cast<double*>(_initialTractions);
constraintVertices = const_cast<int*>(_constraintVertices);
numConstraintVert = _numConstraintVert;
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.hh 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.hh 2011-09-30 23:02:10 UTC (rev 18998)
@@ -67,7 +67,7 @@
static const double _orientation[]; ///< Expected values for fault orientation.
static const double _area[]; ///< Expected values for fault area.
- static const double _forcesInitial[]; ///< Expected values for initial forces.
+ static const double _initialTractions[]; ///< Expected values for initial tractions.
static const double _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const double _slipSlipE[]; ///< Expected values for slip for slip case.
static const double _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/hex8_initialtract.spatialdb
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/hex8_initialtract.spatialdb 2011-09-30 02:42:09 UTC (rev 18997)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/hex8_initialtract.spatialdb 2011-09-30 23:02:10 UTC (rev 18998)
@@ -11,7 +11,7 @@
space-dim = 3
}
}
-0.0 -0.57735027 -0.57735027 1.0000 2.0000 -3.0000
-0.0 0.57735027 -0.57735027 1.1000 2.1000 -3.1000
-0.0 -0.57735027 0.57735027 1.2000 2.2000 -3.2000
-0.0 0.57735027 0.57735027 1.3000 2.3000 -3.3000
+0.0 -1.0 -1.0 1.0000 2.0000 -3.0000
+0.0 1.0 -1.0 1.1000 2.1000 -3.1000
+0.0 -1.0 1.0 1.2000 2.2000 -3.2000
+0.0 1.0 1.0 1.3000 2.3000 -3.3000
More information about the CIG-COMMITS
mailing list