[cig-commits] [commit] knepley/fix-faults-parallel: Ignore clamped edges in various operations. (31e8b9d)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Apr 18 13:33:40 PDT 2014
Repository : ssh://geoshell/pylith
On branch : knepley/fix-faults-parallel
Link : https://github.com/geodynamics/pylith/compare/1cbca9cbd832376cceca629383ad3c8e3db090f0...1169098c7387a0574706ddb12645c08f3401a304
>---------------------------------------------------------------
commit 31e8b9db7cd3784b5a578fef69bc6ae47842f575
Author: Brad Aagaard <baagaard at usgs.gov>
Date: Fri Apr 18 12:44:45 2014 -0700
Ignore clamped edges in various operations.
>---------------------------------------------------------------
31e8b9db7cd3784b5a578fef69bc6ae47842f575
libsrc/pylith/faults/FaultCohesiveLagrange.cc | 84 ++++++++++++++++++++++++---
libsrc/pylith/faults/FaultCohesiveLagrange.hh | 9 +++
2 files changed, 86 insertions(+), 7 deletions(-)
diff --git a/libsrc/pylith/faults/FaultCohesiveLagrange.cc b/libsrc/pylith/faults/FaultCohesiveLagrange.cc
index e5cdd7a..2fb25f2 100644
--- a/libsrc/pylith/faults/FaultCohesiveLagrange.cc
+++ b/libsrc/pylith/faults/FaultCohesiveLagrange.cc
@@ -248,7 +248,10 @@ pylith::faults::FaultCohesiveLagrange::integrateResidual(const topology::Field&
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
- if (e_lagrange < 0) continue;
+ if (e_lagrange < 0) { // Skip clamped edges.
+ continue;
+ } // if
+
// Compute contribution only if Lagrange constraint is local.
PetscInt goff = 0;
err = PetscSectionGetOffset(residualGlobalSection, e_lagrange, &goff);PYLITH_CHECK_ERROR(err);
@@ -389,7 +392,10 @@ pylith::faults::FaultCohesiveLagrange::integrateJacobian(topology::Jacobian* jac
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
- if (e_lagrange < 0) continue;
+ if (e_lagrange < 0) { // Skip clamped edges.
+ continue;
+ } // if
+
// Compute contribution only if Lagrange constraint is local.
PetscInt gloff = 0;
err = PetscSectionGetOffset(solnGlobalSection, e_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
@@ -529,7 +535,10 @@ pylith::faults::FaultCohesiveLagrange::integrateJacobian(topology::Field* jacobi
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
- if (e_lagrange < 0) continue;
+ if (e_lagrange < 0) { // Skip clamped edges
+ continue;
+ } // if
+
// Compute contribution only if Lagrange constraint is local.
PetscInt goff = 0;
err = PetscSectionGetOffset(jacobianGlobalSection, e_lagrange, &goff);PYLITH_CHECK_ERROR(err);
@@ -646,6 +655,10 @@ pylith::faults::FaultCohesiveLagrange::calcPreconditioner(PetscMat* const precon
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
+ if (e_lagrange < 0) { // Skip clamped edges.
+ continue;
+ } // if
+
// Compute contribution only if Lagrange constraint is local.
PetscInt gloff = 0;
err = PetscSectionGetOffset(solnGlobalSection, e_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
@@ -813,6 +826,10 @@ pylith::faults::FaultCohesiveLagrange::adjustSolnLumped(topology::SolutionFields
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
+ if (e_lagrange < 0) { // Skip clamped edges
+ continue;
+ } // if
+
// Set Lagrange multiplier value. Value from preliminary solve is
// bogus due to artificial diagonal entry.
const PetscInt dtloff = dispTIncrVisitor.sectionOffset(e_lagrange);
@@ -1023,6 +1040,12 @@ pylith::faults::FaultCohesiveLagrange::checkConstraints(const topology::Field& s
for(int iVertex = 0; iVertex < numVertices; ++iVertex) {
const int v_negative = _cohesiveVertices[iVertex].negative;
+ const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
+
+ if (e_lagrange < 0) { // Skip clamped edges.
+ continue;
+ } // if
+
assert(spaceDim == solutionVisitor.sectionDof(v_negative));
numConstraints = solutionVisitor.sectionConstraintDof(v_negative);
if (numConstraints > 0) {
@@ -1313,6 +1336,29 @@ pylith::faults::FaultCohesiveLagrange::globalToFault(topology::Field* field,
PYLITH_METHOD_END;
} // faultToGlobal
+
+// ----------------------------------------------------------------------
+// Check to see if given vertex is clamped.
+bool
+pylith::faults::FaultCohesiveLagrange::isClampedVertex(PetscDMLabel clamped,
+ PetscInt vertex)
+{ // _isClampedVertex
+ PYLITH_METHOD_BEGIN;
+
+ bool isClamped = false;
+
+ if (clamped) {
+ PetscInt value = -1;
+ PetscErrorCode err = DMLabelGetValue(clamped, vertex, &value);PYLITH_CHECK_ERROR(err);
+ if (value >= 0) {
+ isClamped = true;
+ } // if
+ } // if
+
+ PYLITH_METHOD_RETURN(isClamped);
+} // _isClampedVertex
+
+
// ----------------------------------------------------------------------
// Calculate orientation at fault vertices.
void
@@ -1353,6 +1399,9 @@ pylith::faults::FaultCohesiveLagrange::_calcOrientation(const PylithScalar upDir
const PetscInt cStart = cellsStratum.begin();
const PetscInt cEnd = cellsStratum.end();
+ PetscDMLabel clamped = NULL;
+ PetscErrorCode err = DMPlexGetLabel(faultDMMesh, "clamped", &clamped);PYLITH_CHECK_ERROR(err);
+
// Containers for orientation information.
// Allocate orientation field.
scalar_array orientationVertex(orientationSize);
@@ -1396,19 +1445,31 @@ pylith::faults::FaultCohesiveLagrange::_calcOrientation(const PylithScalar upDir
// Get orientations at fault cell's vertices.
coordsVisitor.getClosure(&coordsCell, cell);
- PetscErrorCode err = DMPlexGetTransitiveClosure(faultDMMesh, cell, PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
+ err = DMPlexGetTransitiveClosure(faultDMMesh, cell, PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
// Filter out non-vertices
PetscInt numVertices = 0;
+ bool hasClampedVertex = false;
for(PetscInt p = 0; p < closureSize*2; p += 2) {
if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
+ // :KLUDGE: Filter out clamped vertices. Remove this when fault edges are clamped.
+ if (isClampedVertex(clamped, closure[p])) {
+ hasClampedVertex = true;
+ } // if
+
closure[numVertices*2] = closure[p];
closure[numVertices*2+1] = closure[p+1];
++numVertices;
} // if
} // for
+ if (hasClampedVertex) {
+ err = DMPlexRestoreTransitiveClosure(faultDMMesh, cell, PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
+ continue;
+ } // if
for(PetscInt v = 0; v < numVertices; ++v) {
+ const PetscInt v_fault = closure[v*2];
+
// Compute Jacobian and determinant of Jacobian at vertex
cellGeometry.jacobian(&jacobian, &jacobianDet, &coordsCell[0], numBasis, spaceDim, &verticesRef[v*cohesiveDim], cohesiveDim);
@@ -1416,7 +1477,7 @@ pylith::faults::FaultCohesiveLagrange::_calcOrientation(const PylithScalar upDir
cellGeometry.orientation(&orientationVertex, jacobian, jacobianDet, up);
// Update orientation
- const PetscInt ooff = orientationVisitor.sectionOffset(closure[v*2]);
+ const PetscInt ooff = orientationVisitor.sectionOffset(v_fault);
for(PetscInt d = 0; d < orientationSize; ++d) {
orientationArray[ooff+d] += orientationVertex[d];
@@ -1445,6 +1506,11 @@ pylith::faults::FaultCohesiveLagrange::_calcOrientation(const PylithScalar upDir
orientationArray = orientationVisitor.localArray();
int count = 0;
for(PetscInt v = vStart; v < vEnd; ++v, ++count) {
+ // :KLUDGE: Filter out clamped vertices. Remove this when fault edges are clamped.
+ if (isClampedVertex(clamped, v)) {
+ continue;
+ } // if
+
assert(orientationSize == orientationVisitor.sectionDof(v));
const PetscInt ooff = orientationVisitor.sectionOffset(v);
for(PetscInt d = 0; d < orientationSize; ++d) {
@@ -1739,7 +1805,10 @@ pylith::faults::FaultCohesiveLagrange::_calcTractionsChange(topology::Field* tra
const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
- if (e_lagrange < 0) continue;
+ if (e_lagrange < 0) { // skip clamped edges
+ continue;
+ } // if
+
const PetscInt dtoff = dispTVisitor.sectionOffset(e_lagrange);
assert(spaceDim == dispTVisitor.sectionDof(e_lagrange));
@@ -1913,6 +1982,7 @@ pylith::faults::FaultCohesiveLagrange::_getJacobianSubmatrixNP(PetscMat* jacobia
PYLITH_METHOD_END;
} // _getJacobianSubmatrixNP
+
// ----------------------------------------------------------------------
// Get cell field associated with integrator.
const pylith::topology::Field&
@@ -1952,7 +2022,7 @@ pylith::faults::FaultCohesiveLagrange::cellField(const char* name,
msg << "Request for unknown cell field '" << name << "' for fault '" << label() << ".";
throw std::runtime_error(msg.str());
- // Satisfy return values
+ // Satisfy return value
assert(_fields);
const topology::Field& buffer = _fields->get("buffer (vector)");
PYLITH_METHOD_RETURN(buffer);
diff --git a/libsrc/pylith/faults/FaultCohesiveLagrange.hh b/libsrc/pylith/faults/FaultCohesiveLagrange.hh
index 286b662..bfb6723 100644
--- a/libsrc/pylith/faults/FaultCohesiveLagrange.hh
+++ b/libsrc/pylith/faults/FaultCohesiveLagrange.hh
@@ -227,6 +227,15 @@ public :
void globalToFault(topology::Field* field,
const topology::Field& faultOrientation);
+ /** Check to see if given vertex is clamped.
+ *
+ * @param clamped Label in fault mesh associated with clamped vertices.
+ * @param vertex Point in fault mesh.
+ */
+ static
+ bool isClampedVertex(PetscDMLabel clamped,
+ PetscInt vertex);
+
// PROTECTED STRUCTS //////////////////////////////////////////////////
protected :
More information about the CIG-COMMITS
mailing list