[cig-commits] [commit] knepley/upgrade-petsc-interface: Merge branch 'knepley/upgrade-petsc-interface' of github.com:geodynamics/pylith into knepley/upgrade-petsc-interface (0e5e1c1)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Oct 30 16:28:05 PDT 2013
Repository : ssh://geoshell/pylith
On branch : knepley/upgrade-petsc-interface
Link : https://github.com/geodynamics/pylith/compare/ab98cb7666192b03ba35bc19af358487ecd56ce0...0e5e1c14aac776849b92647b93cd74cabed9c592
>---------------------------------------------------------------
commit 0e5e1c14aac776849b92647b93cd74cabed9c592
Merge: 0c6f872 ab98cb7
Author: Brad Aagaard <baagaard at usgs.gov>
Date: Wed Oct 30 16:30:12 2013 -0700
Merge branch 'knepley/upgrade-petsc-interface' of github.com:geodynamics/pylith into knepley/upgrade-petsc-interface
# By Matthew G. Knepley
# Via Matthew G. Knepley
* 'knepley/upgrade-petsc-interface' of github.com:geodynamics/pylith:
Faults: Changing from Lagrange vertices to hybrid edges
Conflicts:
libsrc/pylith/faults/FaultCohesiveLagrange.cc
>---------------------------------------------------------------
0e5e1c14aac776849b92647b93cd74cabed9c592
libsrc/pylith/faults/FaultCohesiveLagrange.cc | 161 +++++++++++++-------------
libsrc/pylith/faults/FaultCohesiveLagrange.hh | 13 ++-
2 files changed, 89 insertions(+), 85 deletions(-)
diff --cc libsrc/pylith/faults/FaultCohesiveLagrange.cc
index 7cc6fed,705b03a..843cafb
--- a/libsrc/pylith/faults/FaultCohesiveLagrange.cc
+++ b/libsrc/pylith/faults/FaultCohesiveLagrange.cc
@@@ -150,12 -148,12 +150,12 @@@ pylith::faults::FaultCohesiveLagrange::
const int numVertices = _cohesiveVertices.size();
for(PetscInt iVertex = 0; iVertex < numVertices; ++iVertex) {
-- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
++ const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
PetscInt dof;
-- err = PetscSectionGetDof(fieldSection, v_lagrange, &dof);PYLITH_CHECK_ERROR(err);assert(spaceDim == dof);
-- err = PetscSectionSetFieldDof(fieldSection, v_lagrange, 0, 0);PYLITH_CHECK_ERROR(err);
-- err = PetscSectionSetFieldDof(fieldSection, v_lagrange, 1, dof);PYLITH_CHECK_ERROR(err);
++ err = PetscSectionGetDof(fieldSection, e_lagrange, &dof);PYLITH_CHECK_ERROR(err);assert(spaceDim == dof);
++ err = PetscSectionSetFieldDof(fieldSection, e_lagrange, 0, 0);PYLITH_CHECK_ERROR(err);
++ err = PetscSectionSetFieldDof(fieldSection, e_lagrange, 1, dof);PYLITH_CHECK_ERROR(err);
} // for
PYLITH_METHOD_END;
@@@ -230,14 -228,14 +230,14 @@@ pylith::faults::FaultCohesiveLagrange::
// Loop over fault vertices
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
-- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
++ const int e_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;
// Compute contribution only if Lagrange constraint is local.
PetscInt goff = 0;
-- err = PetscSectionGetOffset(residualGlobalSection, v_lagrange, &goff);PYLITH_CHECK_ERROR(err);
++ err = PetscSectionGetOffset(residualGlobalSection, e_lagrange, &goff);PYLITH_CHECK_ERROR(err);
if (goff < 0)
continue;
@@@ -261,8 -259,8 +261,8 @@@
const PetscInt dtpoff = dispTVisitor.sectionOffset(v_positive);
assert(spaceDim == dispTVisitor.sectionDof(v_positive));
-- const PetscInt dtloff = dispTVisitor.sectionOffset(v_lagrange);
-- assert(spaceDim == dispTVisitor.sectionDof(v_lagrange));
++ const PetscInt dtloff = dispTVisitor.sectionOffset(e_lagrange);
++ assert(spaceDim == dispTVisitor.sectionDof(e_lagrange));
// Get dispIncr(t->t+dt) at conventional vertices and Lagrange vertex.
const PetscInt dinoff = dispTIncrVisitor.sectionOffset(v_negative);
@@@ -271,8 -269,8 +271,8 @@@
const PetscInt dipoff = dispTIncrVisitor.sectionOffset(v_positive);
assert(spaceDim == dispTIncrVisitor.sectionDof(v_positive));
-- const PetscInt diloff = dispTIncrVisitor.sectionOffset(v_lagrange);
-- assert(spaceDim == dispTIncrVisitor.sectionDof(v_lagrange));
++ const PetscInt diloff = dispTIncrVisitor.sectionOffset(e_lagrange);
++ assert(spaceDim == dispTIncrVisitor.sectionDof(e_lagrange));
// Assemble contributions into field
const PetscInt rnoff = residualVisitor.sectionOffset(v_negative);
@@@ -281,8 -279,8 +281,8 @@@
const PetscInt rpoff = residualVisitor.sectionOffset(v_positive);
assert(spaceDim == residualVisitor.sectionDof(v_positive));
-- const PetscInt rloff = residualVisitor.sectionOffset(v_lagrange);
-- assert(spaceDim == residualVisitor.sectionDof(v_lagrange));
++ const PetscInt rloff = residualVisitor.sectionOffset(e_lagrange);
++ assert(spaceDim == residualVisitor.sectionDof(e_lagrange));
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
@@@ -370,14 -368,14 +370,14 @@@ pylith::faults::FaultCohesiveLagrange::
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
-- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
++ const int e_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;
// Compute contribution only if Lagrange constraint is local.
PetscInt gloff = 0;
-- err = PetscSectionGetOffset(solnGlobalSection, v_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
++ err = PetscSectionGetOffset(solnGlobalSection, e_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
if (gloff < 0)
continue;
@@@ -512,19 -510,19 +512,19 @@@ pylith::faults::FaultCohesiveLagrange::
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
-- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
++ const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
// Compute contribution only if Lagrange constraint is local.
PetscInt goff = 0;
-- err = PetscSectionGetOffset(jacobianGlobalSection, v_lagrange, &goff);PYLITH_CHECK_ERROR(err);
++ err = PetscSectionGetOffset(jacobianGlobalSection, e_lagrange, &goff);PYLITH_CHECK_ERROR(err);
if (goff < 0)
continue;
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(updateEvent);
#endif
-- const PetscInt off = jacobianVisitor.sectionOffset(v_lagrange);
-- assert(spaceDim == jacobianVisitor.sectionDof(v_lagrange));
++ const PetscInt off = jacobianVisitor.sectionOffset(e_lagrange);
++ assert(spaceDim == jacobianVisitor.sectionDof(e_lagrange));
for(PetscInt d = 0; d < spaceDim; ++d) {
jacobianArray[off+d] = 1.0;
@@@ -625,14 -623,14 +625,14 @@@ pylith::faults::FaultCohesiveLagrange::
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0, cV = 0; iVertex < numVertices; ++iVertex) {
-- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
++ const int e_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;
// Compute contribution only if Lagrange constraint is local.
PetscInt gloff = 0;
-- err = PetscSectionGetOffset(solnGlobalSection, v_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
++ err = PetscSectionGetOffset(solnGlobalSection, e_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
if (gloff < 0) {
continue;
} // if
@@@ -697,7 -695,7 +697,7 @@@
INSERT_VALUES);
#if 0 // DEBUGGING
-- std::cout << "1/P_vertex " << *v_lagrange << std::endl;
++ std::cout << "1/P_vertex " << *e_lagrange << std::endl;
for(int iDim = 0; iDim < spaceDim; ++iDim) {
std::cout << " " << precondVertexL[iDim] << std::endl;
} // for
@@@ -792,22 -790,22 +792,22 @@@ pylith::faults::FaultCohesiveLagrange::
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
-- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
++ const int e_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;
// Set Lagrange multiplier value. Value from preliminary solve is
// bogus due to artificial diagonal entry.
-- const PetscInt dtloff = dispTIncrVisitor.sectionOffset(v_lagrange);
-- assert(spaceDim == dispTIncrVisitor.sectionDof(v_lagrange));
++ const PetscInt dtloff = dispTIncrVisitor.sectionOffset(e_lagrange);
++ assert(spaceDim == dispTIncrVisitor.sectionDof(e_lagrange));
for(PetscInt d = 0; d < spaceDim; ++d) {
dispTIncrArray[dtloff+d] = 0.0;
} // for
// Compute contribution only if Lagrange constraint is local.
PetscInt goff;
-- err = PetscSectionGetOffset(jacobianGlobalSection, v_lagrange, &goff);PYLITH_CHECK_ERROR(err);
++ err = PetscSectionGetOffset(jacobianGlobalSection, e_lagrange, &goff);PYLITH_CHECK_ERROR(err);
if (goff < 0) {
continue;
} // if
@@@ -817,8 -815,8 +817,8 @@@
#endif
// Residual at Lagrange vertex.
-- const PetscInt rloff = residualVisitor.sectionOffset(v_lagrange);
-- assert(spaceDim == residualVisitor.sectionDof(v_lagrange));
++ const PetscInt rloff = residualVisitor.sectionOffset(e_lagrange);
++ assert(spaceDim == residualVisitor.sectionDof(e_lagrange));
// Get jacobian at cohesive cell's vertices.
const PetscInt jnoff = jacobianVisitor.sectionOffset(v_negative);
@@@ -845,8 -843,8 +845,8 @@@
const PetscInt dapoff = dispTIncrAdjVisitor.sectionOffset(v_positive);
assert(spaceDim == dispTIncrAdjVisitor.sectionDof(v_positive));
-- const PetscInt daloff = dispTIncrAdjVisitor.sectionOffset(v_lagrange);
-- assert(spaceDim == dispTIncrAdjVisitor.sectionDof(v_lagrange));
++ const PetscInt daloff = dispTIncrAdjVisitor.sectionOffset(e_lagrange);
++ assert(spaceDim == dispTIncrAdjVisitor.sectionDof(e_lagrange));
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
@@@ -888,18 -886,20 +888,21 @@@ pylith::faults::FaultCohesiveLagrange::
PYLITH_METHOD_BEGIN;
assert(_quadrature);
+ PetscErrorCode err = 0;
PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
- topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
- const PetscInt vStart = verticesStratum.begin();
- const PetscInt vEnd = verticesStratum.end();
+ topology::Stratum edgeStratum(dmMesh, topology::Stratum::DEPTH, 1);
+ const PetscInt eEnd = edgeStratum.end();
+ PetscInt eMax;
+
+ err = DMPlexGetHybridBounds(dmMesh, NULL, NULL, &eMax, NULL);PYLITH_CHECK_ERROR(err);
+ // Check for fault groups
PetscBool hasLabel;
- PetscErrorCode err = DMPlexHasLabel(dmMesh, label(), &hasLabel);PYLITH_CHECK_ERROR(err);
+ err = DMPlexHasLabel(dmMesh, label(), &hasLabel);PYLITH_CHECK_ERROR(err);
if (!hasLabel) {
std::ostringstream msg;
- msg << "Mesh missing group of vertices '" << label() << " for boundary condition.";
+ msg << "Mesh missing group of vertices '" << label() << " defining fault.";
throw std::runtime_error(msg.str());
} // if
@@@ -1030,21 -1020,22 +1033,21 @@@ void pylith::faults::FaultCohesiveLagra
assert(_quadrature);
- // :TODO: Update this.
-
- const int numConstraintVert = _quadrature->numBasis();
- const int numCorners = 2 * numConstraintVert; // cohesive cell
- const int numConstraintEdges = _quadrature->numBasis();
++ const int numBasis = _quadrature->numBasis();
+ PetscErrorCode err = 0;
// Get cohesive cells
PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
delete _cohesiveIS; _cohesiveIS = new topology::StratumIS(dmMesh, "material-id", id());assert(_cohesiveIS);
-- const PetscInt* cells = _cohesiveIS->points();
-- const int ncells = _cohesiveIS->size();
++ const PetscInt* cohesiveCells = _cohesiveIS->points();
++ const int numCohesiveCells = _cohesiveIS->size();
-- // Get domain vertices
- topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
- const PetscInt vStart = verticesStratum.begin();
- const PetscInt vEnd = verticesStratum.end();
++ // Get domain edges (Lagrange multiplier DOF)
+ topology::Stratum edgeStratum(dmMesh, topology::Stratum::DEPTH, 1);
+ const PetscInt eStart = edgeStratum.begin();
+ const PetscInt eEnd = edgeStratum.end();
+ PetscInt eMax;
-
+ err = DMPlexGetHybridBounds(dmMesh, NULL, NULL, &eMax, NULL);PYLITH_CHECK_ERROR(err);
// Get vertices and cells in fault mesh.
PetscDM faultDMMesh = _faultMesh->dmMesh();assert(faultDMMesh);
@@@ -1054,7 -1045,7 +1057,7 @@@
topology::Stratum faultCellsStratum(faultDMMesh, topology::Stratum::HEIGHT, 1);
const PetscInt fcStart = faultCellsStratum.begin();
const PetscInt fcEnd = faultCellsStratum.end();
-- assert(ncells == fcEnd-fcStart);
++ assert(numCohesiveCells == fcEnd-fcStart);
topology::SubMeshIS faultMeshIS(*_faultMesh);
const PetscInt numPoints = faultMeshIS.size();
@@@ -1066,54 -1057,56 +1069,54 @@@
_cohesiveVertices.resize(fvEnd-fvStart);
PetscInt index = 0;
- PetscErrorCode err = 0;
-- for(PetscInt c = 0; c < ncells; ++c) {
-- _cohesiveToFault[cells[c]] = c+fcStart;
++ for (PetscInt iCell = 0; iCell < numCohesiveCells; ++iCell) {
++ _cohesiveToFault[cohesiveCells[iCell]] = iCell+fcStart;
// Get oriented closure
PetscInt *closure = NULL;
PetscInt closureSize, q = 0;
-- err = DMPlexGetTransitiveClosure(dmMesh, cells[c], PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
-- for(PetscInt p = 0; p < closureSize*2; p += 2) {
++ err = DMPlexGetTransitiveClosure(dmMesh, cohesiveCells[iCell], PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
++ for (PetscInt p = 0; p < closureSize*2; p += 2) {
const PetscInt point = closure[p];
- if ((point >= vStart) && (point < vEnd)) {
+ if ((point >= eMax) && (point < eEnd)) {
closure[q++] = point;
} // if
} // for
closureSize = q;
- assert(closureSize == numCorners);
- assert(closureSize == numConstraintEdges);
++ assert(closureSize == numBasis);
- for (int iConstraint = 0; iConstraint < numConstraintVert; ++iConstraint) {
- // normal cohesive vertices i and j and constraint vertex k
- const int indexI = iConstraint;
- const int indexJ = iConstraint + numConstraintVert;
- const int indexK = iConstraint + 2 * numConstraintVert;
- for (int iConstraint = 0; iConstraint < numConstraintEdges; ++iConstraint) {
- // normal cohesive vertices i and j and constraint vertex k
- const PetscInt v_lagrange = closure[iConstraint];
- const PetscInt *cone;
- PetscInt coneSize;
++ for (int iConstraint = 0; iConstraint < numBasis; ++iConstraint) {
++ const PetscInt e_lagrange = closure[iConstraint];
- const PetscInt v_lagrange = closure[indexK];
- const PetscInt v_negative = closure[indexI];
- const PetscInt v_positive = closure[indexJ];
- err = DMPlexGetConeSize(dmMesh, v_lagrange, &coneSize);PYLITH_CHECK_ERROR(err);
++ const PetscInt *cone = NULL;
++ PetscInt coneSize;
++ err = DMPlexGetConeSize(dmMesh, e_lagrange, &coneSize);PYLITH_CHECK_ERROR(err);
+ assert(coneSize == 2);
- err = DMPlexGetCone(dmMesh, v_lagrange, &cone);PYLITH_CHECK_ERROR(err);
-
++ err = DMPlexGetCone(dmMesh, e_lagrange, &cone);PYLITH_CHECK_ERROR(err);
+ const PetscInt v_negative = cone[0];
+ const PetscInt v_positive = cone[1];
- PetscInt v_fault;
+ PetscInt v_fault;
err = PetscFindInt(v_negative, numPoints, points, &v_fault);PYLITH_CHECK_ERROR(err);
assert(v_fault >= 0);
-- if (indexMap.end() == indexMap.find(v_lagrange)) {
-- _cohesiveVertices[index].lagrange = v_lagrange;
++ if (indexMap.end() == indexMap.find(e_lagrange)) {
++ _cohesiveVertices[index].lagrange = e_lagrange;
_cohesiveVertices[index].positive = v_positive;
_cohesiveVertices[index].negative = v_negative;
_cohesiveVertices[index].fault = v_fault;
#if 0
- std::cout << "cohesiveVertices[" << index << "]: "
- << "l: " << v_lagrange
+ std::cout << "cohesiveVertices[" << index << "]: "
- << "l: " << v_lagrange
++ << "l: " << e_lagrange
<< ", p: " << v_positive
<< ", n: " << v_negative
<< ", f: " << v_fault
<< std::endl;
#endif
-- indexMap[v_lagrange] = index; // add index to map
++ indexMap[e_lagrange] = index; // add index to map
++index;
} // if
} // for
-- err = DMPlexRestoreTransitiveClosure(dmMesh, cells[c], PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
++ err = DMPlexRestoreTransitiveClosure(dmMesh, cohesiveCells[iCell], PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
} // for
PYLITH_METHOD_END;
@@@ -1660,11 -1651,11 +1663,11 @@@ pylith::faults::FaultCohesiveLagrange::
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
-- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
++ const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
-- const PetscInt dtoff = dispTVisitor.sectionOffset(v_lagrange);
-- assert(spaceDim == dispTVisitor.sectionDof(v_lagrange));
++ const PetscInt dtoff = dispTVisitor.sectionOffset(e_lagrange);
++ assert(spaceDim == dispTVisitor.sectionDof(e_lagrange));
const PetscInt ooff = orientationVisitor.sectionOffset(v_fault);
assert(spaceDim*spaceDim == orientationVisitor.sectionDof(v_fault));
@@@ -1770,11 -1762,11 +1773,11 @@@ pylith::faults::FaultCohesiveLagrange::
const int numVertices = _cohesiveVertices.size();
int numIndicesNP = 0;
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
-- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
++ const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
// Compute contribution only if Lagrange constraint is local.
PetscInt goff = 0;
-- err = PetscSectionGetOffset(solutionGlobalSection, v_lagrange, &goff);PYLITH_CHECK_ERROR(err);
++ err = PetscSectionGetOffset(solutionGlobalSection, e_lagrange, &goff);PYLITH_CHECK_ERROR(err);
if (goff < 0)
continue;
@@@ -1783,13 -1775,13 +1786,13 @@@
int_array indicesNP(numIndicesNP*spaceDim);
for (int iVertex=0, indexNP=0; iVertex < numVertices; ++iVertex) {
-- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
++ const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
// Compute contribution only if Lagrange constraint is local.
PetscInt gloff = 0;
-- err = PetscSectionGetOffset(solutionGlobalSection, v_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
++ err = PetscSectionGetOffset(solutionGlobalSection, e_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
if (gloff < 0)
continue;
diff --cc libsrc/pylith/faults/FaultCohesiveLagrange.hh
index 514c237,514c237..ade1309
--- a/libsrc/pylith/faults/FaultCohesiveLagrange.hh
+++ b/libsrc/pylith/faults/FaultCohesiveLagrange.hh
@@@ -229,14 -229,14 +229,15 @@@ public
// PROTECTED STRUCTS //////////////////////////////////////////////////
protected :
-- /** Data structure to hold relations between vertices in cohesive cell
-- * in mesh of domain and cell of fault mesh.
++ /** Data structure to hold relations between points (vertices and
++ * edges) in cohesive cell in mesh of domain and cell of fault
++ * mesh.
*/
struct CohesiveInfo {
-- int lagrange; ///< Vertex associated with Lagrange multiplier.
-- int positive; ///< Vertex on positive side of the fault.
-- int negative; ///< Vertex on negative side of the fault.
-- int fault; ///< Vertex in fault mesh.
++ int lagrange; ///< Point (edge) associated with Lagrange multiplier.
++ int positive; ///< Point (vertex) on positive side of the fault.
++ int negative; ///< Point (vertex) on negative side of the fault.
++ int fault; ///< Point (vertex) in fault mesh.
};
// PROTECTED METHODS //////////////////////////////////////////////////
More information about the CIG-COMMITS
mailing list