[cig-commits] r21478 - short/3D/PyLith/trunk/libsrc/pylith/faults
brad at geodynamics.org
brad at geodynamics.org
Fri Mar 8 16:12:38 PST 2013
Author: brad
Date: 2013-03-08 16:12:38 -0800 (Fri, 08 Mar 2013)
New Revision: 21478
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
Log:
Code cleanup.
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc 2013-03-08 23:45:29 UTC (rev 21477)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc 2013-03-09 00:12:38 UTC (rev 21478)
@@ -46,11 +46,6 @@
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveImpulses::FaultCohesiveImpulses(void) :
_threshold(1.0e-6),
@@ -196,20 +191,17 @@
PylithScalar scale = 0.0;
int fiberDim = 0;
if (0 == strcasecmp("slip", name)) {
- const topology::Field<topology::SubMesh>& dispRel =
- _fields->get("relative disp");
+ const topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
_allocateBufferVectorField();
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
buffer.copy(dispRel);
buffer.label("slip");
FaultCohesiveLagrange::globalToFault(&buffer, orientation);
return buffer;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
- PetscSection orientationSection = _fields->get("orientation").petscSection();
- Vec orientationVec = _fields->get("orientation").localVector();
- assert(orientationSection);assert(orientationVec);
+ PetscSection orientationSection = _fields->get("orientation").petscSection();assert(orientationSection);
+ PetscVec orientationVec = _fields->get("orientation").localVector();assert(orientationVec);
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
buffer.copy(orientationSection, 0, PETSC_DETERMINE, orientationVec);
@@ -218,8 +210,8 @@
return buffer;
} else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
- PetscSection orientationSection = _fields->get("orientation").petscSection();
- Vec orientationVec = _fields->get("orientation").localVector();
+ PetscSection orientationSection = _fields->get("orientation").petscSection();assert(orientationSection);
+ PetscVec orientationVec = _fields->get("orientation").localVector();assert(orientationVec);
assert(orientationSection);assert(orientationVec);
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
@@ -229,8 +221,8 @@
return buffer;
} else if (0 == strcasecmp("normal_dir", name)) {
- PetscSection orientationSection = _fields->get("orientation").petscSection();
- Vec orientationVec = _fields->get("orientation").localVector();
+ PetscSection orientationSection = _fields->get("orientation").petscSection();assert(orientationSection);
+ PetscVec orientationVec = _fields->get("orientation").localVector();assert(orientationVec);
assert(orientationSection);assert(orientationVec);
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
@@ -297,21 +289,19 @@
amplitude.allocate();
amplitude.scale(lengthScale);
- PetscSection amplitudeSection = amplitude.petscSection();
- Vec amplitudeVec = amplitude.localVector();
- PetscScalar *amplitudeArray;
- assert(amplitudeSection);assert(amplitudeVec);
+ PetscSection amplitudeSection = amplitude.petscSection();assert(amplitudeSection);
+ PetscVec amplitudeVec = amplitude.localVector();assert(amplitudeVec);
+ PetscScalar *amplitudeArray = NULL;
const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
assert(cs);
- DM faultDMMesh = _faultMesh->dmMesh();
- assert(faultDMMesh);
+ PetscDM faultDMMesh = _faultMesh->dmMesh();assert(faultDMMesh);
scalar_array coordsVertex(spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- PetscScalar *coords;
+ PetscSection coordSection = NULL;
+ PetscVec coordVec = NULL;
+ PetscScalar *coordArray = NULL;
err = DMPlexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
@@ -320,23 +310,26 @@
const char* valueNames[1] = { "slip" };
_dbImpulseAmp->queryVals(valueNames, 1);
+ err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
+
std::map<int, int> pointOrder;
int count = 0;
const int numVertices = _cohesiveVertices.size();
- err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
- err = VecGetArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_fault = _cohesiveVertices[iVertex].fault;
- PetscInt cdof, coff;
+ PetscInt cdof, coff;
err = PetscSectionGetDof(coordSection, v_fault, &cdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(coordSection, v_fault, &coff);CHECK_PETSC_ERROR(err);
assert(cdof == spaceDim);
- for(PetscInt d = 0; d < cdof; ++d) coordsVertex[d] = coords[coff+d];
+ for(PetscInt d = 0; d < cdof; ++d) {
+ coordsVertex[d] = coordArray[coff+d];
+ } // for
_normalizer->dimensionalize(&coordsVertex[0], coordsVertex.size(), lengthScale);
- PetscInt adof, aoff;
+ PetscInt adof, aoff;
err = PetscSectionGetDof(amplitudeSection, v_fault, &adof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(amplitudeSection, v_fault, &aoff);CHECK_PETSC_ERROR(err);
assert(adof == fiberDim);
@@ -363,7 +356,7 @@
++count;
} // if
} // for
- err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
// Close properties database
@@ -383,14 +376,13 @@
// Order of impulses is set by processor rank and order of points in
// mesh, using only those points with nonzero amplitudes.
- DM faultDMMesh = _faultMesh->dmMesh();
- PetscErrorCode err;
- assert(faultDMMesh);
+ PetscDM faultDMMesh = _faultMesh->dmMesh();assert(faultDMMesh);
// Gather number of points on each processor.
const int numImpulsesLocal = pointOrder.size();
MPI_Comm comm;
PetscMPIInt commSize, commRank;
+ PetscErrorCode err;
err = PetscObjectGetComm((PetscObject) faultDMMesh, &comm);CHECK_PETSC_ERROR(err);
err = MPI_Comm_size(comm, &commSize);CHECK_PETSC_ERROR(err);
err = MPI_Comm_rank(comm, &commRank);CHECK_PETSC_ERROR(err);
@@ -417,7 +409,7 @@
} // for
} // for
-#if 0 // DEBUGGING
+#if 0 // DEBUGGING :TODO: Update for DM mesh.
const ALE::Obj<RealSection>& amplitudeSection = _fields->get("impulse amplitude").section();
assert(!amplitudeSection.isNull());
int impulse = 0;
@@ -450,21 +442,21 @@
const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
assert(cs);
const int spaceDim = cs->spaceDim();
+
PetscErrorCode err;
- PetscSection amplitudeSection = _fields->get("impulse amplitude").petscSection();
- Vec amplitudeVec = _fields->get("impulse amplitude").localVector();
- PetscScalar *amplitudeArray;
- assert(amplitudeSection);assert(amplitudeVec);
+ PetscSection amplitudeSection = _fields->get("impulse amplitude").petscSection();assert(amplitudeSection);
+ PetscVec amplitudeVec = _fields->get("impulse amplitude").localVector();assert(amplitudeVec);
+ PetscScalar *amplitudeArray = NULL;
- PetscSection dispRelSection = dispRel.petscSection();
- Vec dispRelVec = dispRel.localVector();
- PetscScalar *dispRelArray;
- assert(dispRelSection);assert(dispRelVec);
+ PetscSection dispRelSection = dispRel.petscSection();assert(dispRelSection);
+ PetscVec dispRelVec = dispRel.localVector();assert(dispRelVec);
+ PetscScalar *dispRelArray = NULL;
- const srcs_type::const_iterator& impulseInfo = _impulsePoints.find(impulse);
err = VecGetArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(dispRelVec, &dispRelArray);CHECK_PETSC_ERROR(err);
+
+ const srcs_type::const_iterator& impulseInfo = _impulsePoints.find(impulse);
if (impulseInfo != _impulsePoints.end()) {
const int iVertex = impulseInfo->second.indexCohesive;
const int v_fault = _cohesiveVertices[iVertex].fault;
@@ -472,18 +464,20 @@
// Get amplitude of slip impulse
PetscInt adof, aoff;
-
err = PetscSectionGetDof(amplitudeSection, v_fault, &adof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(amplitudeSection, v_fault, &aoff);CHECK_PETSC_ERROR(err);
assert(adof == 1);
- PetscInt drdof, droff;
+ PetscInt drdof, droff;
err = PetscSectionGetDof(dispRelSection, v_fault, &drdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispRelSection, v_fault, &droff);CHECK_PETSC_ERROR(err);
assert(drdof == spaceDim);
- for(PetscInt d = 0; d < drdof; ++d) {dispRelArray[droff+d] = 0.0;}
- const int indexDOF = impulseInfo->second.indexDOF;
+ for(PetscInt d = 0; d < drdof; ++d) {
+ dispRelArray[droff+d] = 0.0;
+ } // for
+
+ const int indexDOF = impulseInfo->second.indexDOF;
assert(indexDOF >= 0 && indexDOF < spaceDim);
dispRelArray[droff+indexDOF] = amplitudeArray[aoff];
} // if
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveKin.cc 2013-03-08 23:45:29 UTC (rev 21477)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveKin.cc 2013-03-09 00:12:38 UTC (rev 21478)
@@ -47,11 +47,6 @@
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void)
{ // constructor
@@ -99,16 +94,16 @@
pylith::faults::FaultCohesiveKin::initialize(const topology::Mesh& mesh,
const PylithScalar upDir[3])
{ // initialize
- assert(0 != upDir);
- assert(0 != _quadrature);
- assert(0 != _normalizer);
+ assert(upDir);
+ assert(_quadrature);
+ assert(_normalizer);
FaultCohesiveLagrange::initialize(mesh, upDir);
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);
+ assert(src);
src->initialize(*_faultMesh, *_normalizer);
} // for
} // initialize
@@ -122,9 +117,9 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateResidual
- assert(0 != fields);
- assert(0 != _fields);
- assert(0 != _logger);
+ assert(fields);
+ assert(_fields);
+ assert(_logger);
const int setupEvent = _logger->eventId("FaIR setup");
_logger->eventBegin(setupEvent);
@@ -135,7 +130,7 @@
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);
+ assert(src);
if (t >= src->originTime())
src->slip(&dispRel, t);
} // for
@@ -157,10 +152,10 @@
pylith::faults::FaultCohesiveKin::vertexField(const char* name,
const topology::SolutionFields* fields)
{ // vertexField
- assert(0 != _faultMesh);
- assert(0 != _quadrature);
- assert(0 != _normalizer);
- assert(0 != _fields);
+ assert(_faultMesh);
+ assert(_quadrature);
+ assert(_normalizer);
+ assert(_fields);
const int cohesiveDim = _faultMesh->dimension();
const int spaceDim = _quadrature->spaceDim();
@@ -181,9 +176,8 @@
FaultCohesiveLagrange::globalToFault(&buffer, orientation);
return buffer;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
- PetscSection orientationSection = _fields->get("orientation").petscSection();
- Vec orientationVec = _fields->get("orientation").localVector();
- assert(orientationSection);assert(orientationVec);
+ PetscSection orientationSection = _fields->get("orientation").petscSection();assert(orientationSection);
+ PetscVec orientationVec = _fields->get("orientation").localVector();assert(orientationVec);
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
buffer.copy(orientationSection, 0, PETSC_DETERMINE, orientationVec);
@@ -191,8 +185,8 @@
buffer.scale(1.0);
return buffer;
} else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
- PetscSection orientationSection = _fields->get("orientation").petscSection();
- Vec orientationVec = _fields->get("orientation").localVector();
+ PetscSection orientationSection = _fields->get("orientation").petscSection();assert(orientationSection);
+ PetscVec orientationVec = _fields->get("orientation").localVector();assert(orientationVec);
assert(orientationSection);assert(orientationVec);
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
@@ -201,8 +195,8 @@
buffer.scale(1.0);
return buffer;
} else if (0 == strcasecmp("normal_dir", name)) {
- PetscSection orientationSection = _fields->get("orientation").petscSection();
- Vec orientationVec = _fields->get("orientation").localVector();
+ PetscSection orientationSection = _fields->get("orientation").petscSection();assert(orientationSection);
+ PetscVec orientationVec = _fields->get("orientation").localVector();assert(orientationVec);
assert(orientationSection);assert(orientationVec);
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
@@ -245,7 +239,7 @@
return buffer;
} else if (0 == strcasecmp("traction_change", name)) {
- assert(0 != fields);
+ assert(fields);
const topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
@@ -263,7 +257,7 @@
throw std::logic_error("Unknown field in FaultCohesiveKin::vertexField().");
// Satisfy return values
- assert(0 != _fields);
+ assert(_fields);
const topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
return buffer;
} // vertexField
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-03-08 23:45:29 UTC (rev 21477)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-03-09 00:12:38 UTC (rev 21478)
@@ -81,14 +81,14 @@
pylith::faults::FaultCohesiveLagrange::initialize(const topology::Mesh& mesh,
const PylithScalar upDir[3])
{ // initialize
- assert(0 != upDir);
- assert(0 != _quadrature);
- assert(0 != _normalizer);
+ assert(upDir);
+ assert(_quadrature);
+ assert(_normalizer);
_initializeLogger();
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
- assert(0 != cs);
+ assert(cs);
delete _faultMesh;
_faultMesh = new topology::SubMesh();
@@ -114,6 +114,10 @@
logger.stagePop();
+ _quadrature->initializeGeometry();
+
+#if defined(PRECOMPUTE_GEOMETRY)
+#error("Code for PRECOMPUTE_GEOMETRY not implemented.");
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
@@ -122,7 +126,6 @@
const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
_quadrature->initializeGeometry();
-#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->computeGeometry(*_faultMesh, cells);
#endif
@@ -138,52 +141,51 @@
void
pylith::faults::FaultCohesiveLagrange::splitField(topology::Field<topology::Mesh>* field)
{ // splitField
- assert(0 != field);
+ assert(field);
- DM dm = field->dmMesh();
- PetscSection section = field->petscSection();
+ PetscDM dmMesh = field->dmMesh();assert(dmMesh);
+ PetscSection fieldSection = field->petscSection();assert(fieldSection);
const PetscInt spaceDim = field->mesh().dimension();
- PetscInt numFields, numComp, pStart, pEnd;
- PetscErrorCode err;
+ PetscInt numFields, numComp, pStart, pEnd;
- err = PetscSectionGetNumFields(section, &numFields);CHECK_PETSC_ERROR(err);
+ PetscErrorCode err;
+ err = PetscSectionGetNumFields(fieldSection, &numFields);CHECK_PETSC_ERROR(err);
// TODO: Does this make sense?
- if (!numFields) return;
+ if (!numFields)
+ return;
assert(numFields == 2);
- err = PetscSectionGetFieldComponents(section, 0, &numComp);CHECK_PETSC_ERROR(err);
- assert(numComp == spaceDim);
- err = PetscSectionGetFieldComponents(section, 1, &numComp);CHECK_PETSC_ERROR(err);
- assert(numComp == spaceDim);
+ err = PetscSectionGetFieldComponents(fieldSection, 0, &numComp);CHECK_PETSC_ERROR(err);assert(numComp == spaceDim);
+ err = PetscSectionGetFieldComponents(fieldSection, 1, &numComp);CHECK_PETSC_ERROR(err);assert(numComp == spaceDim);
const int numVertices = _cohesiveVertices.size();
for(PetscInt iVertex = 0; iVertex < numVertices; ++iVertex) {
- PetscInt dof;
-
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
- err = PetscSectionGetDof(section, v_lagrange, &dof);CHECK_PETSC_ERROR(err);
+
+ PetscInt dof;
+ err = PetscSectionGetDof(fieldSection, v_lagrange, &dof);CHECK_PETSC_ERROR(err);
assert(spaceDim == dof);
- err = PetscSectionSetFieldDof(section, v_lagrange, 1, dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetFieldDof(fieldSection, v_lagrange, 1, dof);CHECK_PETSC_ERROR(err);
} // for
- err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetChart(fieldSection, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
for(PetscInt p = pStart; p < pEnd; ++p) {
PetscInt dof;
-
- err = PetscSectionGetFieldDof(section, p, 1, &dof);CHECK_PETSC_ERROR(err);
- if (!dof) {err = PetscSectionSetFieldDof(section, p, 0, spaceDim);CHECK_PETSC_ERROR(err);}
+ err = PetscSectionGetFieldDof(fieldSection, p, 1, &dof);CHECK_PETSC_ERROR(err);
+ if (!dof) {
+ err = PetscSectionSetFieldDof(fieldSection, p, 0, spaceDim);CHECK_PETSC_ERROR(err);
+ } // if
}
} // splitField
// ----------------------------------------------------------------------
// Integrate contribution of cohesive cells to residual term.
void
-pylith::faults::FaultCohesiveLagrange::integrateResidual(
- const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::faults::FaultCohesiveLagrange::integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateResidual
- assert(0 != fields);
- assert(0 != _fields);
- assert(0 != _logger);
+ assert(fields);
+ assert(_fields);
+ assert(_logger);
// Cohesive cells with conventional vertices N and P, and constraint
// vertex L make contributions to the assembled residual:
@@ -206,50 +208,46 @@
const int spaceDim = _quadrature->spaceDim();
// Get sections associated with cohesive cells
- DM residualDM = residual.dmMesh();
- PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
- PetscSection residualGlobalSection;
- PetscScalar *residualArray;
+ PetscDM residualDM = residual.dmMesh();assert(residualDM);
+ PetscSection residualSection = residual.petscSection();assert(residualSection);
+ PetscVec residualVec = residual.localVector();assert(residualVec);
+ PetscSection residualGlobalSection = NULL;
+ PetscScalar *residualArray = NULL;
PetscErrorCode err;
- assert(residualSection);assert(residualVec);
err = DMGetDefaultGlobalSection(residualDM, &residualGlobalSection);CHECK_PETSC_ERROR(err);
- PetscSection dispTSection = fields->get("disp(t)").petscSection();
- Vec dispTVec = fields->get("disp(t)").localVector();
- PetscScalar *dispTArray;
- assert(dispTSection);assert(dispTVec);
+ PetscSection dispTSection = fields->get("disp(t)").petscSection();assert(dispTSection);
+ PetscVec dispTVec = fields->get("disp(t)").localVector();assert(dispTVec);
+ PetscScalar *dispTArray = NULL;
- PetscSection dispTIncrSection = fields->get("dispIncr(t->t+dt)").petscSection();
- Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
- PetscScalar *dispTIncrArray;
- assert(dispTIncrSection);assert(dispTIncrVec);
+ PetscSection dispTIncrSection = fields->get("dispIncr(t->t+dt)").petscSection();assert(dispTIncrSection);
+ PetscVec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();assert(dispTIncrVec);
+ PetscScalar *dispTIncrArray = NULL;
- PetscSection dispRelSection = _fields->get("relative disp").petscSection();
- Vec dispRelVec = _fields->get("relative disp").localVector();
- PetscScalar *dispRelArray;
- assert(dispRelSection);assert(dispRelVec);
+ PetscSection dispRelSection = _fields->get("relative disp").petscSection();assert(dispRelSection);
+ PetscVec dispRelVec = _fields->get("relative disp").localVector();assert(dispRelVec);
+ PetscScalar *dispRelArray = NULL;
- PetscSection areaSection = _fields->get("area").petscSection();
- Vec areaVec = _fields->get("area").localVector();
- PetscScalar *areaArray;
- assert(areaSection);assert(areaVec);
+ PetscSection areaSection = _fields->get("area").petscSection();assert(areaSection);
+ PetscVec areaVec = _fields->get("area").localVector();assert(areaVec);
+ PetscScalar *areaArray = NULL;
// Get fault information
- DM dmMesh = fields->mesh().dmMesh();
+ PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
#endif
- // Loop over fault vertices
- const int numVertices = _cohesiveVertices.size();
err = VecGetArray(dispRelVec, &dispRelArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(areaVec, &areaArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
+
+ // 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 v_fault = _cohesiveVertices[iVertex].fault;
@@ -267,80 +265,69 @@
// Get relative dislplacement at fault vertex.
PetscInt drdof, droff;
-
err = PetscSectionGetDof(dispRelSection, v_fault, &drdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispRelSection, v_fault, &droff);CHECK_PETSC_ERROR(err);
assert(spaceDim == drdof);
// Get area associated with fault vertex.
PetscInt adof, aoff;
-
err = PetscSectionGetDof(areaSection, v_fault, &adof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(areaSection, v_fault, &aoff);CHECK_PETSC_ERROR(err);
assert(1 == adof);
+ const PylithScalar area = areaArray[aoff];
// Get disp(t) at conventional vertices and Lagrange vertex.
PetscInt dtndof, dtnoff;
-
err = PetscSectionGetDof(dispTSection, v_negative, &dtndof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispTSection, v_negative, &dtnoff);CHECK_PETSC_ERROR(err);
assert(spaceDim == dtndof);
- PetscInt dtpdof, dtpoff;
+ PetscInt dtpdof, dtpoff;
err = PetscSectionGetDof(dispTSection, v_positive, &dtpdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispTSection, v_positive, &dtpoff);CHECK_PETSC_ERROR(err);
assert(spaceDim == dtpdof);
- PetscInt dtldof, dtloff;
+ PetscInt dtldof, dtloff;
err = PetscSectionGetDof(dispTSection, v_lagrange, &dtldof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispTSection, v_lagrange, &dtloff);CHECK_PETSC_ERROR(err);
assert(spaceDim == dtldof);
// Get dispIncr(t->t+dt) at conventional vertices and Lagrange vertex.
PetscInt dindof, dinoff;
-
err = PetscSectionGetDof(dispTIncrSection, v_negative, &dindof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispTIncrSection, v_negative, &dinoff);CHECK_PETSC_ERROR(err);
assert(spaceDim == dindof);
- PetscInt dipdof, dipoff;
+ PetscInt dipdof, dipoff;
err = PetscSectionGetDof(dispTIncrSection, v_positive, &dipdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispTIncrSection, v_positive, &dipoff);CHECK_PETSC_ERROR(err);
assert(spaceDim == dipdof);
- PetscInt dildof, diloff;
+ PetscInt dildof, diloff;
err = PetscSectionGetDof(dispTIncrSection, v_lagrange, &dildof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispTIncrSection, v_lagrange, &diloff);CHECK_PETSC_ERROR(err);
assert(spaceDim == dildof);
-#if defined(DETAILED_EVENT_LOGGING)
- _logger->eventEnd(restrictEvent);
- _logger->eventBegin(computeEvent);
-#endif
-
- // Compute current estimate of displacement at time t+dt using solution increment.
- for(PetscInt d = 0; d < spaceDim; ++d) {
- } // for
-
-#if defined(DETAILED_EVENT_LOGGING)
- _logger->eventEnd(computeEvent);
- _logger->eventBegin(updateEvent);
-#endif
-
// Assemble contributions into field
PetscInt rndof, rnoff, rpdof, rpoff, rldof, rloff;
-
err = PetscSectionGetDof(residualSection, v_negative, &rndof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(residualSection, v_negative, &rnoff);CHECK_PETSC_ERROR(err);
err = PetscSectionGetDof(residualSection, v_positive, &rpdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(residualSection, v_positive, &rpoff);CHECK_PETSC_ERROR(err);
err = PetscSectionGetDof(residualSection, v_lagrange, &rldof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(residualSection, v_lagrange, &rloff);CHECK_PETSC_ERROR(err);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(restrictEvent);
+ _logger->eventBegin(updateEvent);
+#endif
+
assert(spaceDim == rndof);assert(spaceDim == rpdof);assert(spaceDim == rldof);
for(PetscInt d = 0; d < spaceDim; ++d) {
- residualArray[rnoff+d] += areaArray[aoff] * (dispTArray[dtloff+d] + dispTIncrArray[diloff+d]);
- residualArray[rpoff+d] += -areaArray[aoff] * (dispTArray[dtloff+d] + dispTIncrArray[diloff+d]);
- residualArray[rloff+d] += -areaArray[aoff] * (dispTArray[dtpoff+d] + dispTIncrArray[dipoff+d] - dispTArray[dtnoff+d] - dispTIncrArray[dinoff+d] - dispRelArray[droff+d]);
+ const PylithScalar residualN = area * (dispTArray[dtloff+d] + dispTIncrArray[diloff+d]);
+ residualArray[rnoff+d] += +residualN;
+ residualArray[rpoff+d] += -residualN;
+ residualArray[rloff+d] += -area * (dispTArray[dtpoff+d] + dispTIncrArray[dipoff+d] - dispTArray[dtnoff+d] - dispTIncrArray[dinoff+d] - dispRelArray[droff+d]);
}
#if defined(DETAILED_EVENT_LOGGING)
@@ -362,15 +349,14 @@
// ----------------------------------------------------------------------
// Compute Jacobian matrix (A) associated with operator.
void
-pylith::faults::FaultCohesiveLagrange::integrateJacobian(
- topology::Jacobian* jacobian,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::faults::FaultCohesiveLagrange::integrateJacobian(topology::Jacobian* jacobian,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateJacobian
- assert(0 != jacobian);
- assert(0 != fields);
- assert(0 != _fields);
- assert(0 != _logger);
+ assert(jacobian);
+ assert(fields);
+ assert(_fields);
+ assert(_logger);
const int setupEvent = _logger->eventId("FaIJ setup");
const int geometryEvent = _logger->eventId("FaIJ geometry");
@@ -388,13 +374,13 @@
// Get sections
PetscSection areaSection = _fields->get("area").petscSection();
- Vec areaVec = _fields->get("area").localVector();
+ PetscVec areaVec = _fields->get("area").localVector();
PetscScalar *areaArray;
assert(areaSection);assert(areaVec);
- DM solnDM = fields->solution().dmMesh();
+ PetscDM solnDM = fields->solution().dmMesh();
PetscSection solnSection = fields->solution().petscSection();
- Vec solnVec = fields->solution().localVector();
+ PetscVec solnVec = fields->solution().localVector();
PetscSection solnGlobalSection;
PetscScalar *solnArray;
PetscErrorCode err;
@@ -402,7 +388,7 @@
err = DMGetDefaultGlobalSection(solnDM, &solnGlobalSection);CHECK_PETSC_ERROR(err);
// Get fault information
- DM dmMesh = fields->mesh().dmMesh();
+ PetscDM dmMesh = fields->mesh().dmMesh();
assert(dmMesh);
// Allocate vectors for vertex values
@@ -416,7 +402,7 @@
// Get sparse matrix
const PetscMat jacobianMatrix = jacobian->matrix();
- assert(0 != jacobianMatrix);
+ assert(jacobianMatrix);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -537,10 +523,10 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateJacobian
- assert(0 != jacobian);
- assert(0 != fields);
- assert(0 != _fields);
- assert(0 != _logger);
+ assert(jacobian);
+ assert(fields);
+ assert(_fields);
+ assert(_logger);
const int setupEvent = _logger->eventId("FaIJ setup");
const int geometryEvent = _logger->eventId("FaIJ geometry");
@@ -556,9 +542,9 @@
// multipliers as part of the solve.
const PetscInt spaceDim = _quadrature->spaceDim();
- DM jacobianDM = jacobian->dmMesh();
+ PetscDM jacobianDM = jacobian->dmMesh();
PetscSection jacobianSection = jacobian->petscSection();
- Vec jacobianVec = jacobian->localVector();
+ PetscVec jacobianVec = jacobian->localVector();
PetscSection jacobianGlobalSection;
PetscScalar *jacobianArray;
PetscErrorCode err;
@@ -661,18 +647,18 @@
indicesRel[i] = i;
// Get sections
- DM areaDM = _fields->get("area").dmMesh();
+ PetscDM areaDM = _fields->get("area").dmMesh();
PetscSection areaSection = _fields->get("area").petscSection();
- Vec areaVec = _fields->get("area").localVector();
+ PetscVec areaVec = _fields->get("area").localVector();
PetscSection areaGlobalSection;
PetscScalar *areaArray;
PetscErrorCode err;
assert(areaSection);assert(areaVec);
err = DMGetDefaultGlobalSection(areaDM, &areaGlobalSection);CHECK_PETSC_ERROR(err);
- DM solutionDM = fields->solution().dmMesh();
+ PetscDM solutionDM = fields->solution().dmMesh();
PetscSection solutionSection = fields->solution().petscSection();
- Vec solutionVec = fields->solution().localVector();
+ PetscVec solutionVec = fields->solution().localVector();
PetscSection solutionGlobalSection;
PetscScalar *solutionArray;
assert(solutionSection);assert(solutionVec);
@@ -791,8 +777,8 @@
const topology::Field<
topology::Mesh>& jacobian)
{ // adjustSolnLumped
- assert(0 != fields);
- assert(0 != _quadrature);
+ assert(fields);
+ assert(_quadrature);
// Cohesive cells with conventional vertices N and P, and constraint
// vertex L require 2 adjustments to the solution:
@@ -826,13 +812,13 @@
// Get section information
PetscSection areaSection = _fields->get("area").petscSection();
- Vec areaVec = _fields->get("area").localVector();
+ PetscVec areaVec = _fields->get("area").localVector();
PetscScalar *areaArray;
assert(areaSection);assert(areaVec);
- DM jacobianDM = jacobian.dmMesh();
+ PetscDM jacobianDM = jacobian.dmMesh();
PetscSection jacobianSection = jacobian.petscSection();
- Vec jacobianVec = jacobian.localVector();
+ PetscVec jacobianVec = jacobian.localVector();
PetscSection jacobianGlobalSection;
PetscScalar *jacobianArray;
PetscErrorCode err;
@@ -840,17 +826,17 @@
err = DMGetDefaultGlobalSection(jacobianDM, &jacobianGlobalSection);CHECK_PETSC_ERROR(err);
PetscSection residualSection = fields->get("residual").petscSection();
- Vec residualVec = fields->get("residual").localVector();
+ PetscVec residualVec = fields->get("residual").localVector();
PetscScalar *residualArray;
assert(residualSection);assert(residualVec);
PetscSection dispTIncrSection = fields->get("dispIncr(t->t+dt)").petscSection();
- Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ PetscVec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
PetscScalar *dispTIncrArray;
assert(dispTIncrSection);assert(dispTIncrVec);
PetscSection dispTIncrAdjSection = fields->get("dispIncr adjust").petscSection();
- Vec dispTIncrAdjVec = fields->get("dispIncr adjust").localVector();
+ PetscVec dispTIncrAdjVec = fields->get("dispIncr adjust").localVector();
PetscScalar *dispTIncrAdjArray;
assert(dispTIncrAdjSection);assert(dispTIncrAdjVec);
@@ -963,9 +949,9 @@
void
pylith::faults::FaultCohesiveLagrange::verifyConfiguration(const topology::Mesh& mesh) const
{ // verifyConfiguration
- assert(0 != _quadrature);
+ assert(_quadrature);
- DM dmMesh = mesh.dmMesh();
+ PetscDM dmMesh = mesh.dmMesh();
PetscBool hasLabel;
PetscInt vStart, vEnd;
PetscErrorCode err;
@@ -1095,10 +1081,10 @@
// Initialize auxiliary cohesive cell information.
void pylith::faults::FaultCohesiveLagrange::_initializeCohesiveInfo(const topology::Mesh& mesh)
{ // _initializeCohesiveInfo
- assert(0 != _quadrature);
+ assert(_quadrature);
// Get cohesive cells
- DM dmMesh = mesh.dmMesh();
+ PetscDM dmMesh = mesh.dmMesh();
IS cellIS;
const PetscInt *cells;
PetscInt numCells, vStart, vEnd;
@@ -1111,7 +1097,7 @@
const int numConstraintVert = _quadrature->numBasis();
const int numCorners = 3 * numConstraintVert; // cohesive cell
- DM faultDMMesh = _faultMesh->dmMesh();
+ PetscDM faultDMMesh = _faultMesh->dmMesh();
IS subpointIS;
const PetscInt *points;
PetscInt numPoints, fcStart, fcEnd, fvStart, fvEnd;
@@ -1192,7 +1178,7 @@
{ // initializeLogger
delete _logger;
_logger = new utils::EventLogger;
- assert(0 != _logger);
+ assert(_logger);
_logger->className("FaultCohesiveLagrange");
_logger->initialize();
@@ -1238,16 +1224,16 @@
// Get sections.
PetscSection fieldSection = field->petscSection();
- Vec fieldVec = field->localVector();
+ PetscVec fieldVec = field->localVector();
PetscScalar *fieldArray;
assert(fieldSection);assert(fieldVec);
PetscSection orientationSection = faultOrientation.petscSection();
- Vec orientationVec = faultOrientation.localVector();
+ PetscVec orientationVec = faultOrientation.localVector();
PetscScalar *orientationArray;
assert(orientationSection);assert(orientationVec);
- DM dmMesh = field->mesh().dmMesh();
+ PetscDM dmMesh = field->mesh().dmMesh();
PetscInt vStart, vEnd;
PetscErrorCode err;
@@ -1301,16 +1287,16 @@
// Get sections.
PetscSection fieldSection = field->petscSection();
- Vec fieldVec = field->localVector();
+ PetscVec fieldVec = field->localVector();
PetscScalar *fieldArray;
assert(fieldSection);assert(fieldVec);
PetscSection orientationSection = faultOrientation.petscSection();
- Vec orientationVec = faultOrientation.localVector();
+ PetscVec orientationVec = faultOrientation.localVector();
PetscScalar *orientationArray;
assert(orientationSection);assert(orientationVec);
- DM dmMesh = field->mesh().dmMesh();
+ PetscDM dmMesh = field->mesh().dmMesh();
PetscInt vStart, vEnd;
PetscErrorCode err;
@@ -1352,9 +1338,9 @@
void
pylith::faults::FaultCohesiveLagrange::_calcOrientation(const PylithScalar upDir[3])
{ // _calcOrientation
- assert(0 != upDir);
- assert(0 != _faultMesh);
- assert(0 != _fields);
+ assert(upDir);
+ assert(_faultMesh);
+ assert(_fields);
scalar_array up(3);
for (int i=0; i < 3; ++i)
@@ -1362,7 +1348,7 @@
// Get vertices in fault mesh.
MPI_Comm comm;
- DM faultDMMesh = _faultMesh->dmMesh();
+ PetscDM faultDMMesh = _faultMesh->dmMesh();
PetscInt vStart, vEnd, cStart, cEnd;
PetscErrorCode err;
@@ -1404,7 +1390,7 @@
orientation.allocate();
orientation.zero();
PetscSection orientationSection = orientation.petscSection();
- Vec orientationVec = orientation.localVector();
+ PetscVec orientationVec = orientation.localVector();
PetscScalar *orientationArray;
logger.stagePop();
@@ -1413,7 +1399,7 @@
// Get section containing coordinates of vertices
PetscSection coordSection;
- Vec coordVec;
+ PetscVec coordVec;
err = DMPlexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
@@ -1664,7 +1650,7 @@
scalar_array areaCell(numBasis);
// Get vertices in fault mesh.
- DM faultDMMesh = _faultMesh->dmMesh();
+ PetscDM faultDMMesh = _faultMesh->dmMesh();
PetscInt vStart, vEnd, cStart, cEnd;
PetscErrorCode err;
@@ -1686,7 +1672,7 @@
area.scale(pow(lengthScale, (spaceDim-1)));
area.zero();
PetscSection areaSection = area.petscSection();
- Vec areaVec = area.localVector();
+ PetscVec areaVec = area.localVector();
PetscScalar *areaArray;
assert(areaSection);assert(areaVec);
@@ -1694,7 +1680,7 @@
scalar_array coordinatesCell(numBasis*spaceDim);
PetscSection coordSection;
- Vec coordVec;
+ PetscVec coordVec;
err = DMPlexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
@@ -1747,10 +1733,10 @@
topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& dispT)
{ // _calcTractionsChange
- assert(0 != tractions);
- assert(0 != _faultMesh);
- assert(0 != _fields);
- assert(0 != _normalizer);
+ assert(tractions);
+ assert(_faultMesh);
+ assert(_fields);
+ assert(_normalizer);
tractions->label("traction_change");
tractions->scale(_normalizer->pressureScale());
@@ -1760,13 +1746,13 @@
// Get sections.
PetscSection dispTSection = dispT.petscSection();
- Vec dispTVec = dispT.localVector();
+ PetscVec dispTVec = dispT.localVector();
PetscScalar *dispTArray;
PetscErrorCode err;
assert(dispTSection);assert(dispTVec);
PetscSection orientationSection = _fields->get("orientation").petscSection();
- Vec orientationVec = _fields->get("orientation").localVector();
+ PetscVec orientationVec = _fields->get("orientation").localVector();
PetscScalar *orientationArray;
assert(orientationSection);assert(orientationVec);
@@ -1781,7 +1767,7 @@
tractionsSection = tractions->petscSection();
logger.stagePop();
} // if
- Vec tractionsVec = tractions->localVector();
+ PetscVec tractionsVec = tractions->localVector();
PetscScalar *tractionsArray;
assert(tractionsSection);assert(tractionsVec);
tractions->zero();
@@ -1829,7 +1815,7 @@
void
pylith::faults::FaultCohesiveLagrange::_allocateBufferVectorField(void)
{ // _allocateBufferVectorField
- assert(0 != _fields);
+ assert(_fields);
if (_fields->hasField("buffer (vector)"))
return;
@@ -1838,7 +1824,7 @@
// Create vector field; use same shape/chart as relative
// displacement field.
- assert(0 != _faultMesh);
+ assert(_faultMesh);
_fields->add("buffer (vector)", "buffer");
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
const topology::Field<topology::SubMesh>& dispRel =
@@ -1855,7 +1841,7 @@
void
pylith::faults::FaultCohesiveLagrange::_allocateBufferScalarField(void)
{ // _allocateBufferScalarField
- assert(0 != _fields);
+ assert(_fields);
if (_fields->hasField("buffer (scalar)"))
return;
@@ -1863,7 +1849,7 @@
logger.stagePush("OutputFields");
// Create vector field; use same shape/chart as area field.
- assert(0 != _faultMesh);
+ assert(_faultMesh);
_fields->add("buffer (scalar)", "buffer");
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (scalar)");
buffer.newSection(topology::FieldBase::VERTICES_FIELD, 1);
@@ -1890,9 +1876,9 @@
assert(indicesMatToSubmat);
// Get global order
- DM solutionDM = fields.solution().dmMesh();
+ PetscDM solutionDM = fields.solution().dmMesh();
PetscSection solutionSection = fields.solution().petscSection();
- Vec solutionVec = fields.solution().localVector();
+ PetscVec solutionVec = fields.solution().localVector();
PetscSection solutionGlobalSection;
PetscScalar *solutionArray;
PetscErrorCode err;
@@ -1901,10 +1887,10 @@
// Get Jacobian matrix
const PetscMat jacobianMatrix = jacobian.matrix();
- assert(0 != jacobianMatrix);
+ assert(jacobianMatrix);
const spatialdata::geocoords::CoordSys* cs = fields.mesh().coordsys();
- assert(0 != cs);
+ assert(cs);
const int spaceDim = cs->spaceDim();
const int numVertices = _cohesiveVertices.size();
@@ -1978,7 +1964,7 @@
{ // cellField
if (0 == strcasecmp("partition", name)) {
- DM faultDMMesh = _faultMesh->dmMesh();
+ PetscDM faultDMMesh = _faultMesh->dmMesh();
PetscInt cStart, cEnd;
PetscErrorCode err;
@@ -1993,7 +1979,7 @@
topology::Field<topology::SubMesh>& partition = _fields->get("partition");
partition.allocate();
PetscSection partitionSection = partition.petscSection();
- Vec partitionVec = partition.localVector();
+ PetscVec partitionVec = partition.localVector();
PetscScalar *partitionArray;
assert(partitionSection);assert(partitionVec);
@@ -2025,7 +2011,7 @@
throw std::runtime_error(msg.str());
// Satisfy return values
- assert(0 != _fields);
+ assert(_fields);
const topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
return buffer;
} // cellField
More information about the CIG-COMMITS
mailing list