[cig-commits] r21669 - in short/3D/PyLith/trunk: libsrc/pylith/faults unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 27 17:19:35 PDT 2013
Author: brad
Date: 2013-03-27 17:19:35 -0700 (Wed, 27 Mar 2013)
New Revision: 21669
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc
Log:
Code cleanup. Updated to use visitors.
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2013-03-27 22:43:22 UTC (rev 21668)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2013-03-28 00:19:35 UTC (rev 21669)
@@ -31,6 +31,7 @@
#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
#include "pylith/friction/FrictionModel.hh" // USES FrictionModel
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
#include "pylith/problems/SolverLinear.hh" // USES SolverLinear
@@ -165,7 +166,6 @@
velRel.vectorFieldType(topology::FieldBase::VECTOR);
velRel.scale(_normalizer->lengthScale() / _normalizer->timeScale());
-
PYLITH_METHOD_END;
} // initialize
@@ -202,50 +202,42 @@
// Get cell geometry information that doesn't depend on cell
const int spaceDim = _quadrature->spaceDim();
- PetscErrorCode err = 0;
-
// Get sections associated with cohesive cells
PetscDM residualDM = residual.dmMesh();assert(residualDM);
- PetscSection residualSection = residual.petscSection();assert(residualSection);
- PetscVec residualVec = residual.localVector();assert(residualVec);
PetscSection residualGlobalSection = NULL;
- PetscScalar *residualArray = NULL;
- err = DMGetDefaultGlobalSection(residualDM, &residualGlobalSection);CHECK_PETSC_ERROR(err);
- err = VecGetArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
+ PetscErrorCode err = DMGetDefaultGlobalSection(residualDM, &residualGlobalSection);CHECK_PETSC_ERROR(err);assert(residualGlobalSection);
- PetscSection dispTSection = fields->get("disp(t)").petscSection();assert(dispTSection);
- PetscVec dispTVec = fields->get("disp(t)").localVector();assert(dispTVec);
- PetscScalar *dispTArray = NULL;
- err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
+ topology::VecVisitorMesh residualVisitor(residual);
+ PetscScalar* residualArray = residualVisitor.localArray();
- 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;
- err = VecGetArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ topology::VecVisitorMesh dispTVisitor(dispT);
+ const PetscScalar* dispTArray = dispTVisitor.localArray();
+ topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
+ topology::VecVisitorMesh dispTIncrVisitor(dispTIncr);
+ const PetscScalar* dispTIncrArray = dispTIncrVisitor.localArray();
+
scalar_array tractPerturbVertex(spaceDim);
- PetscSection tractionsSection = NULL;
- PetscVec tractionsVec = NULL;
+ topology::VecVisitorMesh* tractionsVisitor = 0;
PetscScalar *tractionsArray = NULL;
if (_tractPerturbation) {
_tractPerturbation->calculate(t);
- const topology::Fields<topology::Field<topology::SubMesh> >* params = _tractPerturbation->parameterFields();
- assert(params);
- tractionsSection = params->get("value").petscSection();assert(tractionsSection);
- tractionsVec = params->get("value").localVector();assert(tractionsVec);
+ const topology::Fields<topology::Field<topology::SubMesh> >* params = _tractPerturbation->parameterFields();assert(params);
+ const topology::Field<topology::SubMesh>& tractions = params->get("value");
+
+ tractionsVisitor = new topology::VecVisitorMesh(tractions);
+ tractionsArray = tractionsVisitor->localArray();
} // if
- err = VecGetArray(tractionsVec, &tractionsArray);CHECK_PETSC_ERROR(err);
- PetscSection areaSection = _fields->get("area").petscSection();assert(areaSection);
- PetscVec areaVec = _fields->get("area").localVector();assert(areaVec);
- PetscScalar *areaArray = NULL;
- err = VecGetArray(areaVec, &areaArray);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::SubMesh>& area = _fields->get("area");
+ topology::VecVisitorMesh areaVisitor(area);
+ const PetscScalar* areaArray = areaVisitor.localArray();
- PetscSection orientationSection = _fields->get("orientation").petscSection();assert(orientationSection);
- PetscVec orientationVec = _fields->get("orientation").localVector();assert(orientationVec);
- PetscScalar *orientationArray = NULL;
- err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
+ topology::VecVisitorMesh orientationVisitor(area);
+ const PetscScalar* orientationArray = orientationVisitor.localArray();
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -259,11 +251,12 @@
const int v_fault = _cohesiveVertices[iVertex].fault;
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
- PetscInt goff;
// Compute contribution only if Lagrange constraint is local.
+ PetscInt goff = 0;
err = PetscSectionGetOffset(residualGlobalSection, v_lagrange, &goff);CHECK_PETSC_ERROR(err);
- if (goff < 0) continue;
+ if (goff < 0)
+ continue;
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(restrictEvent);
@@ -271,61 +264,42 @@
// Get prescribed traction perturbation at fault vertex.
if (_tractPerturbation) {
- PetscInt vdof, voff;
- err = PetscSectionGetDof(tractionsSection, v_fault, &vdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(tractionsSection, v_fault, &voff);CHECK_PETSC_ERROR(err);
- assert(vdof == spaceDim);
-
+ const PetscInt toff = tractionsVisitor->sectionOffset(v_fault);
+ assert(spaceDim == tractionsVisitor->sectionDof(v_fault));
for(PetscInt d = 0; d < spaceDim; ++d) {
- tractPerturbVertex[d] = tractionsArray[voff+d];
+ tractPerturbVertex[d] = tractionsArray[toff+d];
} // for
} else {
tractPerturbVertex = 0.0;
} // if/else
// Get orientation associated with fault vertex.
- PetscInt odof, ooff;
- err = PetscSectionGetDof(orientationSection, v_fault, &odof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(orientationSection, v_fault, &ooff);CHECK_PETSC_ERROR(err);
- assert(spaceDim*spaceDim == odof);
+ const PetscInt ooff = orientationVisitor.sectionOffset(v_fault);
+ assert(spaceDim*spaceDim == orientationVisitor.sectionDof(v_fault));
// 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 PetscInt aoff = areaVisitor.sectionOffset(v_fault);
+ assert(1 == areaVisitor.sectionDof(v_fault));
// 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;
+ const PetscInt dtnoff = dispTVisitor.sectionOffset(v_negative);
+ assert(spaceDim == dispTVisitor.sectionDof(v_negative));
- err = PetscSectionGetDof(dispTSection, v_positive, &dtpdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(dispTSection, v_positive, &dtpoff);CHECK_PETSC_ERROR(err);
- assert(spaceDim == dtpdof);
+ const PetscInt dtpoff = dispTVisitor.sectionOffset(v_positive);
+ assert(spaceDim == dispTVisitor.sectionDof(v_positive));
- 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);
+ const PetscInt dtloff = dispTVisitor.sectionOffset(v_lagrange);
+ assert(spaceDim == dispTVisitor.sectionDof(v_lagrange));
// 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);
+ const PetscInt dinoff = dispTIncrVisitor.sectionOffset(v_negative);
+ assert(spaceDim == dispTIncrVisitor.sectionDof(v_negative));
- 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);
+ const PetscInt dipoff = dispTIncrVisitor.sectionOffset(v_positive);
+ assert(spaceDim == dispTIncrVisitor.sectionDof(v_positive));
- 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);
+ const PetscInt diloff = dispTIncrVisitor.sectionOffset(v_lagrange);
+ assert(spaceDim == dispTIncrVisitor.sectionDof(v_lagrange));
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
@@ -341,7 +315,6 @@
tractionNormal += orientationArray[ooff+indexN*spaceDim+d] * (dispTArray[dtloff+d] + dispTIncrArray[diloff+d]);
} // for
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
@@ -349,12 +322,12 @@
if (slipNormal < _zeroTolerance || !_openFreeSurf) {
// if no opening or flag indicates to still impose initial tractions when fault is open.
// Assemble contributions into field
- PetscInt rndof, rnoff, rpdof, rpoff;
- 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);
- assert(spaceDim == rndof);assert(spaceDim == rpdof);
+ const PetscInt rnoff = residualVisitor.sectionOffset(v_negative);
+ assert(spaceDim == residualVisitor.sectionDof(v_negative));
+
+ const PetscInt rpoff = residualVisitor.sectionOffset(v_positive);
+ assert(spaceDim == residualVisitor.sectionDof(v_positive));
+
// Initial (external) tractions oppose (internal) tractions associated with Lagrange multiplier.
for(PetscInt d = 0; d < spaceDim; ++d) {
residualArray[rnoff+d] += areaArray[aoff] * (dispTArray[dtloff+d] + dispTIncrArray[diloff+d] - tractPerturbVertex[d]);
@@ -376,12 +349,7 @@
#endif
} // for
PetscLogFlops(numVertices*spaceDim*8);
- err = VecRestoreArray(areaVec, &areaArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(tractionsVec, &tractionsArray);CHECK_PETSC_ERROR(err);
+ delete tractionsVisitor; tractionsVisitor = 0;
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc 2013-03-27 22:43:22 UTC (rev 21668)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc 2013-03-28 00:19:35 UTC (rev 21669)
@@ -30,6 +30,8 @@
#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
@@ -187,7 +189,7 @@
// Get vertex field associated with integrator.
const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::FaultCohesiveImpulses::vertexField(const char* name,
- const topology::SolutionFields* fields)
+ const topology::SolutionFields* fields)
{ // vertexField
PYLITH_METHOD_BEGIN;
@@ -292,8 +294,8 @@
assert(_normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
- const int spaceDim = _quadrature->spaceDim();
- PetscErrorCode err;
+ const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();assert(cs);
+ const int spaceDim = cs->spaceDim();
// Create section to hold amplitudes of impulses.
_fields->add("impulse amplitude", "impulse_amplitude");
@@ -304,50 +306,50 @@
amplitude.allocate();
amplitude.scale(lengthScale);
- PetscSection amplitudeSection = amplitude.petscSection();assert(amplitudeSection);
- PetscVec amplitudeVec = amplitude.localVector();assert(amplitudeVec);
- PetscScalar *amplitudeArray = NULL;
+ topology::VecVisitorMesh amplitudeVisitor(amplitude);
+ PetscScalar *amplitudeArray = amplitudeVisitor.localArray();
- const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
- assert(cs);
-
+ PetscErrorCode err;
PetscDM faultDMMesh = _faultMesh->dmMesh();assert(faultDMMesh);
+ PetscIS vertexNumbering = NULL;
+ err = DMPlexGetVertexNumbering(faultDMMesh, &vertexNumbering);CHECK_PETSC_ERROR(err);
+ const PetscInt* points = NULL;
+ PetscInt npoints = 0;
+ err = ISGetIndices(vertexNumbering, &points);CHECK_PETSC_ERROR(err);
+ err = ISGetLocalSize(vertexNumbering, &npoints);CHECK_PETSC_ERROR(err);
scalar_array coordsVertex(spaceDim);
- 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);
+ topology::CoordsVisitor coordsVisitor(faultDMMesh);
+ PetscScalar *coordsArray = coordsVisitor.localArray();
assert(_dbImpulseAmp);
_dbImpulseAmp->open();
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();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_fault = _cohesiveVertices[iVertex].fault;
- 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);
+ std::cerr << ":TODO: MATT Update FaultCohesiveImpulses::_setupImpulses() for PETSc DM." << std::endl;
+#if 0
+ // Only create impulses on local vertices
+ if (!globalNumbering->isLocal(v_fault)) {
+ continue;
+ } // if
+#endif
- for(PetscInt d = 0; d < cdof; ++d) {
- coordsVertex[d] = coordArray[coff+d];
+ const PetscInt coff = coordsVisitor.sectionOffset(v_fault);
+ assert(spaceDim == coordsVisitor.sectionDof(v_fault));
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coordsVertex[d] = coordsArray[coff+d];
} // for
_normalizer->dimensionalize(&coordsVertex[0], coordsVertex.size(), lengthScale);
- 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);
+ const PetscInt aoff = amplitudeVisitor.sectionOffset(v_fault);
+ assert(fiberDim == amplitudeVisitor.sectionDof(v_fault));
amplitudeArray[aoff] = 0.0;
int err = _dbImpulseAmp->query(&litudeArray[aoff], 1, &coordsVertex[0], coordsVertex.size(), cs);
@@ -371,8 +373,6 @@
++count;
} // if
} // for
- err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
// Close properties database
_dbImpulseAmp->close();
@@ -399,10 +399,9 @@
// Gather number of points on each processor.
const int numImpulsesLocal = pointOrder.size();
- MPI_Comm comm;
+ MPI_Comm comm = _faultMesh->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);
int_array numImpulsesAll(commSize);
@@ -428,22 +427,6 @@
} // for
} // for
-#if 0 // DEBUGGING :TODO: Update for DM mesh.
- const ALE::Obj<RealSection>& amplitudeSection = _fields->get("impulse amplitude").section();
- assert(!amplitudeSection.isNull());
- int impulse = 0;
- for (int irank=0; irank < commSize; ++irank) {
- MPI_Barrier(comm);
- if (commRank == irank) {
- for (int i=0; i < _impulsePoints.size(); ++i, ++impulse) {
- const ImpulseInfoStruct& info = _impulsePoints[impulse];
- const PylithScalar* amplitudeVertex = amplitudeSection->restrictPoint(_cohesiveVertices[info.indexCohesive].fault);
- std::cout << "["<<irank<<"]: " << impulse << " -> (" << info.indexCohesive << "," << info.indexDOF << "), v_fault: " << _cohesiveVertices[info.indexCohesive].fault << ", amplitude: " << amplitudeVertex[0] << std::endl;
- } // for
- } // if
- } // for
-#endif
-
PYLITH_METHOD_END;
} // _setupImpulseOrder
@@ -465,19 +448,13 @@
const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();assert(cs);
const int spaceDim = cs->spaceDim();
- PetscErrorCode err;
+ topology::Field<topology::SubMesh>& amplitude = _fields->get("impulse amplitude");
+ topology::VecVisitorMesh amplitudeVisitor(amplitude);
+ const PetscScalar* amplitudeArray = amplitudeVisitor.localArray();
- 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();assert(dispRelSection);
- PetscVec dispRelVec = dispRel.localVector();assert(dispRelVec);
- PetscScalar *dispRelArray = NULL;
+ topology::VecVisitorMesh dispRelVisitor(dispRel);
+ PetscScalar* dispRelArray = dispRelVisitor.localArray();
- 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;
@@ -485,17 +462,12 @@
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
// 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);
+ const PetscInt aoff = amplitudeVisitor.sectionOffset(v_fault);
+ assert(1 == amplitudeVisitor.sectionDof(v_fault));
- 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) {
+ const PetscInt droff = dispRelVisitor.sectionOffset(v_fault);
+ assert(spaceDim == dispRelVisitor.sectionDof(v_fault));
+ for(PetscInt d = 0; d < spaceDim; ++d) {
dispRelArray[droff+d] = 0.0;
} // for
@@ -503,8 +475,6 @@
assert(indexDOF >= 0 && indexDOF < spaceDim);
dispRelArray[droff+indexDOF] = amplitudeArray[aoff];
} // if
- err = VecRestoreArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(dispRelVec, &dispRelArray);CHECK_PETSC_ERROR(err);
#if 0 // DEBUGGING
std::cout << "impulse: " << impulse << std::endl;
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-03-27 22:43:22 UTC (rev 21668)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-03-28 00:19:35 UTC (rev 21669)
@@ -148,7 +148,7 @@
err = PetscSectionGetNumFields(fieldSection, &numFields);CHECK_PETSC_ERROR(err);
// TODO: Does this make sense?
if (!numFields)
- return;
+ PYLITH_METHOD_END;
assert(numFields == 2);
err = PetscSectionGetFieldComponents(fieldSection, 0, &numComp);CHECK_PETSC_ERROR(err);assert(numComp == spaceDim);
err = PetscSectionGetFieldComponents(fieldSection, 1, &numComp);CHECK_PETSC_ERROR(err);assert(numComp == spaceDim);
@@ -1705,7 +1705,7 @@
assert(_fields);
if (_fields->hasField("buffer (vector)"))
- return;
+ PYLITH_METHOD_END;
// Create vector field; use same shape/chart as relative
// displacement field.
@@ -1730,7 +1730,7 @@
assert(_fields);
if (_fields->hasField("buffer (scalar)"))
- return;
+ PYLITH_METHOD_END;
// Create vector field; use same shape/chart as area field.
assert(_faultMesh);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc 2013-03-27 22:43:22 UTC (rev 21668)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc 2013-03-28 00:19:35 UTC (rev 21669)
@@ -30,6 +30,8 @@
#include <stdlib.h> // USES abort()
+//#define MALLOC_DUMP
+
int
main(int argc,
char* argv[])
@@ -39,7 +41,9 @@
try {
// Initialize PETSc
PetscErrorCode err = PetscInitialize(&argc, &argv, NULL, NULL);CHKERRQ(err);
+#if defined(MALLOC_DUMP)
err = PetscOptionsSetValue("-malloc_dump", "");CHKERRQ(err);
+#endif
// Initialize Python (to eliminate need to initialize when
// parsing units in spatial databases).
@@ -76,6 +80,10 @@
abort();
} // catch
+#if !defined(MALLOC_DUMP)
+ std::cout << "WARNING -malloc dump is OFF\n" << std::endl;
+#endif
+
return (result.wasSuccessful() ? 0 : 1);
} // main
More information about the CIG-COMMITS
mailing list