[cig-commits] r21635 - short/3D/PyLith/trunk/libsrc/pylith/faults
brad at geodynamics.org
brad at geodynamics.org
Mon Mar 25 16:12:17 PDT 2013
Author: brad
Date: 2013-03-25 16:12:17 -0700 (Mon, 25 Mar 2013)
New Revision: 21635
Removed:
short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.hh
short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.icc
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/ConstRateSlipFn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/EqKinSrc.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.hh
short/3D/PyLith/trunk/libsrc/pylith/faults/LiuCosSlipFn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/StepSlipFn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/TimeHistorySlipFn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/faultsfwd.hh
Log:
Code cleanup. Update to use visitors.
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -238,8 +238,6 @@
const PetscInt vStart = depthStratum.begin();
const PetscInt vEnd = depthStratum.end();
- PetscErrorCode err;
-
// Get sections
const topology::Field<topology::SubMesh>& finalSlip = _parameters->get("final slip");
topology::VecVisitorMesh finalSlipVisitor(finalSlip);
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/ConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/ConstRateSlipFn.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/ConstRateSlipFn.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -23,6 +23,9 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/Stratum.hh" // USES Stratum
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -70,20 +73,18 @@
assert(_dbSlipRate);
assert(_dbSlipTime);
- const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
- assert(cs);
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();assert(cs);
const int spaceDim = cs->spaceDim();
const PylithScalar lengthScale = normalizer.lengthScale();
const PylithScalar timeScale = normalizer.timeScale();
- const PylithScalar velocityScale =
- normalizer.lengthScale() / normalizer.timeScale();
+ const PylithScalar velocityScale = normalizer.lengthScale() / normalizer.timeScale();
// Get vertices in fault mesh
PetscDM dmMesh = faultMesh.dmMesh();assert(dmMesh);
- PetscInt vStart, vEnd;
- PetscErrorCode err;
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ topology::Stratum depthStratum(dmMesh, topology::Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
assert(_parameters);
@@ -93,9 +94,8 @@
slipRate.allocate();
slipRate.scale(velocityScale);
slipRate.vectorFieldType(topology::FieldBase::VECTOR);
- PetscSection slipRateSection = slipRate.petscSection();assert(slipRateSection);
- PetscVec slipRateVec = slipRate.localVector();assert(slipRateVec);
- PetscScalar *slipRateArray = NULL;
+ topology::VecVisitorMesh slipRateVisitor(slipRate);
+ PetscScalar* slipRateArray = slipRateVisitor.localArray();
_parameters->add("slip time", "slip_time");
topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
@@ -103,9 +103,8 @@
slipTime.allocate();
slipTime.scale(timeScale);
slipTime.vectorFieldType(topology::FieldBase::SCALAR);
- PetscSection slipTimeSection = slipTime.petscSection();assert(slipTimeSection);
- PetscVec slipTimeVec = slipTime.localVector();assert(slipTimeVec);
- PetscScalar *slipTimeArray = NULL;
+ topology::VecVisitorMesh slipTimeVisitor(slipTime);
+ PetscScalar* slipTimeArray = slipTimeVisitor.localArray();
// Open databases and set query values
_dbSlipRate->open();
@@ -139,31 +138,22 @@
// Get coordinates of vertices
scalar_array vCoordsGlobal(spaceDim);
- PetscSection coordSection = NULL;
- PetscVec coordVec = NULL;
- PetscScalar *coordArray = NULL;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);assert(coordSection);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar* coordsArray = coordsVisitor.localArray();
_slipRateVertex.resize(spaceDim);
- err = VecGetArray(slipRateVec, &slipRateArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt srdof, sroff, stdof, stoff, cdof, coff;
- err = PetscSectionGetDof(slipRateSection, v, &srdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipRateSection, v, &sroff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(coordSection, v, &cdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(coordSection, v, &coff);CHECK_PETSC_ERROR(err);
- assert(cdof == spaceDim);
-
- for(PetscInt d = 0; d < cdof; ++d) {
- vCoordsGlobal[d] = coordArray[coff+d];
+ // Dimensionalize coordinates
+ const PetscInt coff = coordsVisitor.sectionOffset(v);
+ assert(spaceDim == coordsVisitor.sectionDof(v));
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ vCoordsGlobal[d] = coordsArray[coff+d];
} // for
normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(), lengthScale);
+ // Slip Rate
+ const PetscInt sroff = slipRateVisitor.sectionOffset(v);
+ assert(spaceDim == slipRateVisitor.sectionDof(v));
int err = _dbSlipRate->query(&_slipRateVertex[0], _slipRateVertex.size(), &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -175,6 +165,9 @@
} // if
normalizer.nondimensionalize(&_slipRateVertex[0], _slipRateVertex.size(), velocityScale);
+ // Slip time
+ const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+ assert(1 == slipTimeVisitor.sectionDof(v));
err = _dbSlipTime->query(&_slipTimeVertex, 1, &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -188,14 +181,11 @@
// add origin time to rupture time
_slipTimeVertex += originTime;
- for(PetscInt d = 0; d < srdof; ++d) {
+ for(PetscInt d = 0; d < spaceDim; ++d) {
slipRateArray[sroff+d] = _slipRateVertex[d];
} // for
slipTimeArray[stoff] = _slipTimeVertex;
} // for
- err = VecRestoreArray(slipRateVec, &slipRateArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
// Close databases
_dbSlipRate->close();
@@ -212,50 +202,40 @@
assert(_parameters);
// Get vertices in fault mesh
- PetscDM dmMesh = slip->mesh().dmMesh();assert(dmMesh);
- PetscInt vStart, vEnd;
- PetscErrorCode err;
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ PetscDM dmMesh = _parameters->mesh().dmMesh();assert(dmMesh);
+ topology::Stratum depthStratum(dmMesh, topology::Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
// Get sections
const topology::Field<topology::SubMesh>& slipRate = _parameters->get("slip rate");
- PetscSection slipRateSection = slipRate.petscSection();assert(slipRateSection);
- PetscVec slipRateVec = slipRate.localVector();assert(slipRateVec);
- PetscScalar *slipRateArray = NULL;
- err = VecGetArray(slipRateVec, &slipRateArray);CHECK_PETSC_ERROR(err);
+ topology::VecVisitorMesh slipRateVisitor(slipRate);
+ const PetscScalar* slipRateArray = slipRateVisitor.localArray();
const topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
- PetscSection slipTimeSection = slipTime.petscSection();assert(slipTimeSection);
- PetscVec slipTimeVec = slipTime.localVector();assert(slipTimeVec);
- PetscScalar *slipTimeArray = NULL;
- err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
+ topology::VecVisitorMesh slipTimeVisitor(slipTime);
+ const PetscScalar* slipTimeArray = slipTimeVisitor.localArray();
- PetscSection slipSection = slip->petscSection();assert(slipSection);
- PetscVec slipVec = slip->localVector();assert(slipVec);
- PetscScalar *slipArray = NULL;
- err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+ topology::VecVisitorMesh slipVisitor(*slip);
+ PetscScalar* slipArray = slipVisitor.localArray();
+ const int spaceDim = _slipRateVertex.size();
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt srdof, sroff, stdof, stoff, sdof, soff;
- err = PetscSectionGetDof(slipRateSection, v, &srdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipRateSection, v, &sroff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipSection, v, &sdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipSection, v, &soff);CHECK_PETSC_ERROR(err);
- assert(srdof == sdof);
- assert(stdof == 1);
+ const PetscInt sroff = slipRateVisitor.sectionOffset(v);
+ const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+ const PetscInt soff = slipVisitor.sectionOffset(v);
+ assert(spaceDim == slipRateVisitor.sectionDof(v));
+ assert(1 == slipTimeVisitor.sectionDof(v));
+ assert(spaceDim == slipVisitor.sectionDof(v));
+
const PylithScalar relTime = t - slipTimeArray[stoff];
if (relTime > 0.0) {
- for(PetscInt d = 0; d < sdof; ++d) {
+ for(PetscInt d = 0; d < spaceDim; ++d) {
slipArray[soff+d] += slipRateArray[sroff+d] * relTime; // Convert slip rate to slip
} // for
} // if
} // for
- err = VecRestoreArray(slipRateVec, &slipRateArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
PetscLogFlops((vEnd-vStart) * (1 + _slipRateVertex.size()));
} // slip
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/EqKinSrc.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/EqKinSrc.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/EqKinSrc.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -76,9 +76,8 @@
// ----------------------------------------------------------------------
// Initialize slip time function.
void
-pylith::faults::EqKinSrc::initialize(
- const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer)
+pylith::faults::EqKinSrc::initialize(const topology::SubMesh& faultMesh,
+ const spatialdata::units::Nondimensional& normalizer)
{ // initialize
// :TODO: Normalize slip time in Python?
normalizer.nondimensionalize(&_originTime, 1, normalizer.timeScale());
@@ -89,9 +88,8 @@
// ----------------------------------------------------------------------
// Get slip on fault surface at time t.
void
-pylith::faults::EqKinSrc::slip(
- topology::Field<topology::SubMesh>* const slipField,
- const PylithScalar t)
+pylith::faults::EqKinSrc::slip(topology::Field<topology::SubMesh>* const slipField,
+ const PylithScalar t)
{ // slip
assert(_slipfn);
_slipfn->slip(slipField, t);
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -73,14 +73,14 @@
if (!_useFaultMesh) {
// Get group of vertices associated with fault
- DM dmMesh = mesh.dmMesh();
- PetscBool has;
+ PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
+ PetscBool hasLabel;
PetscErrorCode err;
assert(dmMesh);
assert(std::string("") != label());
- err = DMPlexHasLabel(dmMesh, label(), &has);CHECK_PETSC_ERROR(err);
- if (!has) {
+ err = DMPlexHasLabel(dmMesh, label(), &hasLabel);CHECK_PETSC_ERROR(err);
+ if (!hasLabel) {
std::ostringstream msg;
msg << "Mesh missing group of vertices '" << label()
<< "' for fault interface condition.";
@@ -107,13 +107,14 @@
assert(mesh);
assert(std::string("") != label());
+ std::cerr << ":TODO: MATT Update FaultCohesive::adjustTopology for PETSc DM." << std::endl;
+
try {
topology::SubMesh faultMesh;
ALE::Obj<SieveFlexMesh> faultBoundary;
// Get group of vertices associated with fault
- DM dmMesh = mesh->dmMesh();
- assert(dmMesh);
+ PetscDM dmMesh = mesh->dmMesh();assert(dmMesh);
if (!_useFaultMesh) {
PetscDMLabel groupField;
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -43,16 +43,10 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
-#include <petscblaslapack.h> // USES svd and dgemm
-
//#define PRECOMPUTE_GEOMETRY
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveLagrange::FaultCohesiveLagrange(void)
{ // constructor
@@ -79,7 +73,7 @@
// Initialize fault. Determine orientation and setup boundary
void
pylith::faults::FaultCohesiveLagrange::initialize(const topology::Mesh& mesh,
- const PylithScalar upDir[3])
+ const PylithScalar upDir[3])
{ // initialize
assert(upDir);
assert(_quadrature);
@@ -946,12 +940,11 @@
{ // verifyConfiguration
assert(_quadrature);
- PetscDM dmMesh = mesh.dmMesh();
- PetscBool hasLabel;
- PetscInt vStart, vEnd;
+ PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
+ PetscBool hasLabel;
+ PetscInt vStart, vEnd;
PetscErrorCode err;
- assert(dmMesh);
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
err = DMPlexHasLabel(dmMesh, label(), &hasLabel);CHECK_PETSC_ERROR(err);
if (!hasLabel) {
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -20,19 +20,15 @@
#include "FaultCohesiveTract.hh" // implementation of object methods
#include "CohesiveTopology.hh" // USES CohesiveTopology
-#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/topology/Stratum.hh" // USES StratumIS
+
#include <cassert> // USES assert()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::SubMesh::DomainSieveMesh SieveMesh;
-
-// ----------------------------------------------------------------------
// Default constructor.
-pylith::faults::FaultCohesiveTract::FaultCohesiveTract(void) :
- _dbInitial(0)
+pylith::faults::FaultCohesiveTract::FaultCohesiveTract(void)
{ // constructor
_useLagrangeConstraints = false;
} // constructor
@@ -51,17 +47,9 @@
{ // deallocate
FaultCohesive::deallocate();
- _dbInitial = 0; // :TODO: Use shared pointer
} // deallocate
// ----------------------------------------------------------------------
-// Sets the spatial database for the inital tractions
-void pylith::faults::FaultCohesiveTract::dbInitial(spatialdata::spatialdb::SpatialDB* dbs)
-{ // dbInitial
- _dbInitial = dbs;
-} // dbInitial
-
-// ----------------------------------------------------------------------
// Initialize fault. Determine orientation and setup boundary
void
pylith::faults::FaultCohesiveTract::initialize(const topology::Mesh& mesh,
@@ -82,23 +70,14 @@
// Initialize quadrature geometry.
_quadrature->initializeGeometry();
- // Compute orientation at quadrature points in fault mesh.
- _calcOrientation(upDir);
-
- // Get initial tractions using a spatial database.
- _getInitialTractions();
-
- // Setup fault constitutive model.
- _initConstitutiveModel();
} // initialize
// ----------------------------------------------------------------------
// Integrate contribution of cohesive cells to residual term.
void
-pylith::faults::FaultCohesiveTract::integrateResidual(
- const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::faults::FaultCohesiveTract::integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateResidual
throw std::logic_error("FaultCohesiveTract::integrateResidual() not implemented.");
} // integrateResidual
@@ -106,10 +85,9 @@
// ----------------------------------------------------------------------
// Compute Jacobian matrix (A) associated with operator.
void
-pylith::faults::FaultCohesiveTract::integrateJacobian(
- topology::Jacobian* jacobian,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::faults::FaultCohesiveTract::integrateJacobian(topology::Jacobian* jacobian,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateJacobian
throw std::logic_error("FaultCohesiveTract::integrateJacobian() not implemented.");
@@ -119,15 +97,14 @@
// ----------------------------------------------------------------------
// Verify configuration is acceptable.
void
-pylith::faults::FaultCohesiveTract::verifyConfiguration(
- const topology::Mesh& mesh) const
+pylith::faults::FaultCohesiveTract::verifyConfiguration(const topology::Mesh& mesh) const
{ // verifyConfiguration
- assert(0 != _quadrature);
+ assert(_quadrature);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
-
- if (!sieveMesh->hasIntSection(label())) {
+ const PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
+ PetscBool hasLabel = PETSC_FALSE;
+ PetscErrorCode err = DMPlexHasLabel(dmMesh, label(), &hasLabel);CHECK_PETSC_ERROR(err);
+ if (!hasLabel) {
std::ostringstream msg;
msg << "Mesh missing group of vertices '" << label()
<< " for boundary condition.";
@@ -147,20 +124,18 @@
} // if
const int numCorners = _quadrature->refGeometry().numCorners();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", id());
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- const int cellNumCorners = sieveMesh->getNumCellCorners(*c_iter);
- if (2*numCorners != cellNumCorners) {
+ topology::StratumIS cohesiveIS(dmMesh, "material-id", id());
+ const PetscInt* cells = cohesiveIS.points();
+ const PetscInt ncells = cohesiveIS.size();
+
+ PetscInt coneSize = 0;
+ for (PetscInt i=0; i < ncells; ++i) {
+ err = DMPlexGetConeSize(dmMesh, cells[i], &coneSize);CHECK_PETSC_ERROR(err);
+ if (2*numCorners != coneSize) {
std::ostringstream msg;
msg << "Number of vertices in reference cell (" << numCorners
- << ") is not compatible with number of vertices (" << cellNumCorners
- << ") in cohesive cell " << *c_iter << " for fault '"
+ << ") is not compatible with number of vertices (" << coneSize
+ << ") in cohesive cell " << cells[i] << " for fault '"
<< label() << "'.";
throw std::runtime_error(msg.str());
} // if
@@ -170,9 +145,8 @@
// ----------------------------------------------------------------------
// Get vertex field associated with integrator.
const pylith::topology::Field<pylith::topology::SubMesh>&
-pylith::faults::FaultCohesiveTract::vertexField(
- const char* name,
- const topology::SolutionFields* fields)
+pylith::faults::FaultCohesiveTract::vertexField(const char* name,
+ const topology::SolutionFields* fields)
{ // vertexField
throw std::logic_error("FaultCohesiveTract::vertexField() not implemented.");
} // vertexField
@@ -180,254 +154,11 @@
// ----------------------------------------------------------------------
// Get cell field associated with integrator.
const pylith::topology::Field<pylith::topology::SubMesh>&
-pylith::faults::FaultCohesiveTract::cellField(
- const char* name,
- const topology::SolutionFields* fields)
+pylith::faults::FaultCohesiveTract::cellField(const char* name,
+ const topology::SolutionFields* fields)
{ // cellField
throw std::logic_error("FaultCohesiveTract::cellField() not implemented.");
} // cellField
-// ----------------------------------------------------------------------
-// Calculate orientation at fault vertices.
-void
-pylith::faults::FaultCohesiveTract::_calcOrientation(const PylithScalar upDir[3])
-{ // _calcOrientation
- assert(0 != _fields);
- assert(0 != _quadrature);
- scalar_array up(3);
- for (int i=0; i < 3; ++i)
- up[i] = upDir[i];
-
- // Get 'fault' cells.
- DM faultDMMesh = _faultMesh->dmMesh();
- PetscInt cStart, cEnd;
- PetscErrorCode err;
-
- assert(faultDMMesh);
- err = DMPlexGetHeightStratum(faultDMMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
-
- // Quadrature related values.
- const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
- const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
- const int numBasis = _quadrature->numBasis();
- const int numQuadPts = _quadrature->numQuadPts();
- const int spaceDim = cellGeometry.spaceDim();
- scalar_array quadPtRef(cellDim);
- const scalar_array& quadPtsRef = _quadrature->quadPtsRef();
-
- // Containers for orientation information
- const int orientationSize = spaceDim * spaceDim;
- const int fiberDim = numQuadPts * orientationSize;
- const int jacobianSize = spaceDim * cellDim;
- scalar_array jacobian(jacobianSize);
- PylithScalar jacobianDet = 0;
- scalar_array orientationQuadPt(orientationSize);
- scalar_array orientationCell(fiberDim);
-
- // Get sections.
- scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- err = DMPlexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
-
- // :TODO: Use spaces to create subsections like in FaultCohesiveKin.
- _fields->add("orientation", "orientation",
- topology::FieldBase::CELLS_FIELD, fiberDim);
- topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
- orientation.allocate();
- PetscSection orientationSection = orientation.petscSection();
- Vec orientationVec = orientation.localVector();
- PetscScalar *orientationArray;
- assert(orientationSection);assert(orientationVec);
-
- // Loop over cells in fault mesh and compute orientations.
- err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
- for(PetscInt c = cStart; c < cEnd; ++c) {
- PetscScalar *coords = PETSC_NULL;
- PetscInt coordsSize;
-
- err = DMPlexVecGetClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- err = DMPlexVecRestoreClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- _quadrature->computeGeometry(coordinatesCell, c);
-
- // Reset orientation to zero.
- orientationCell = 0.0;
-
- // Compute orientation at each quadrature point of current cell.
- for (int iQuad=0, iRef=0, iSpace=0;
- iQuad < numQuadPts;
- ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
- // Reset orientation at quad pt to zero.
- orientationQuadPt = 0.0;
-
- // Compute Jacobian and determinant at quadrature point, then get
- // orientation.
- memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(PylithScalar));
- cellGeometry.jacobian(&jacobian, &jacobianDet,
- coordinatesCell, quadPtRef);
- cellGeometry.orientation(&orientationQuadPt, jacobian, jacobianDet, up);
- assert(jacobianDet > 0.0);
- orientationQuadPt /= jacobianDet;
-
- memcpy(&orientationCell[iQuad*orientationSize],
- &orientationQuadPt[0], orientationSize*sizeof(PylithScalar));
- } // for
- PetscInt dof, off;
-
- err = PetscSectionGetDof(orientationSection, c, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(orientationSection, c, &off);CHECK_PETSC_ERROR(err);
- for(PetscInt d = 0; d < dof; ++d) {
- orientationArray[off+d] = orientationCell[d];
- }
- } // for
- err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
-
- // debugging
- orientation.view("FAULT ORIENTATION");
-} // _calcOrientation
-
-// ----------------------------------------------------------------------
-void
-pylith::faults::FaultCohesiveTract::_getInitialTractions(void)
-{ // _getInitialTractions
- assert(0 != _normalizer);
- assert(0 != _quadrature);
-
- const PylithScalar pressureScale = _normalizer->pressureScale();
- const PylithScalar lengthScale = _normalizer->lengthScale();
-
- const int spaceDim = _quadrature->spaceDim();
- const int numQuadPts = _quadrature->numQuadPts();
-
- if (0 != _dbInitial) { // Setup initial values, if provided.
- // Create section to hold initial tractions.
- _fields->add("initial traction", "initial_traction");
- topology::Field<topology::SubMesh>& traction = _fields->get("initial traction");
- traction.newSection(topology::FieldBase::CELLS_FIELD, numQuadPts*spaceDim);
- traction.allocate();
- traction.scale(pressureScale);
- traction.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- PetscSection tractionSection = traction.petscSection();
- Vec tractionVec = traction.localVector();
- PetscScalar *tractionArray;
- assert(tractionSection);assert(tractionVec);
-
- _dbInitial->open();
- switch (spaceDim)
- { // switch
- case 1 : {
- const char* valueNames[] = {"traction-normal"};
- _dbInitial->queryVals(valueNames, 1);
- break;
- } // case 1
- case 2 : {
- const char* valueNames[] = {"traction-shear", "traction-normal"};
- _dbInitial->queryVals(valueNames, 2);
- break;
- } // case 2
- case 3 : {
- const char* valueNames[] = {"traction-shear-leftlateral",
- "traction-shear-updip",
- "traction-normal"};
- _dbInitial->queryVals(valueNames, 3);
- break;
- } // case 3
- default :
- std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
- assert(0);
- throw std::logic_error("Bad spatial dimension in Neumann.");
- } // switch
-
- // Get 'fault' cells.
- DM faultDMMesh = _faultMesh->dmMesh();
- PetscInt cStart, cEnd;
- PetscErrorCode err;
-
- assert(faultDMMesh);
- err = DMPlexGetHeightStratum(faultDMMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
-
- const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
- const int numBasis = _quadrature->numBasis();
- const int numQuadPts = _quadrature->numQuadPts();
- const int spaceDim = _quadrature->spaceDim();
-
- // Containers for database query results and quadrature coordinates in
- // reference geometry.
- scalar_array tractionCell(numQuadPts*spaceDim);
- scalar_array quadPtsGlobal(numQuadPts*spaceDim);
-
- // Get sections.
- scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- err = DMPlexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
-
- const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
-
- // Compute quadrature information
-
- // Loop over cells in boundary mesh and perform queries.
- err = VecGetArray(tractionVec, &tractionArray);CHECK_PETSC_ERROR(err);
- for (PetscInt c = cStart; c < cEnd; ++c) {
- PetscScalar *coords = PETSC_NULL;
- PetscInt coordsSize;
-
- err = DMPlexVecGetClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- err = DMPlexVecRestoreClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- _quadrature->computeGeometry(coordinatesCell, c);
-
- const scalar_array& quadPtsNondim = _quadrature->quadPts();
- quadPtsGlobal = quadPtsNondim;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
-
- tractionCell = 0.0;
- for (int iQuad=0, iSpace=0; iQuad < numQuadPts; ++iQuad, iSpace+=spaceDim) {
- const int err = _dbInitial->query(&tractionCell[iQuad*spaceDim], spaceDim,
- &quadPtsGlobal[iSpace], spaceDim, cs);
- if (err) {
- std::ostringstream msg;
- msg << "Could not find initial tractions at (";
- for (int i=0; i < spaceDim; ++i)
- msg << " " << quadPtsGlobal[i+iSpace];
- msg << ") for dynamic fault interface " << label() << "\n"
- << "using spatial database " << _dbInitial->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- } // for
- _normalizer->nondimensionalize(&tractionCell[0], tractionCell.size(),
- pressureScale);
-
- // Update section
- PetscInt dof, off;
-
- err = PetscSectionGetDof(tractionSection, c, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(tractionSection, c, &off);CHECK_PETSC_ERROR(err);
- assert(tractionCell.size() == dof);
- for(PetscInt d = 0; d < dof; ++d) {
- tractionArray[off+d] = tractionCell[d];
- }
- } // for
- err = VecRestoreArray(tractionVec, &tractionArray);CHECK_PETSC_ERROR(err);
-
- _dbInitial->close();
-
- // debugging
- traction.view("INITIAL TRACTIONS");
- } // if
-} // _getInitialTractions
-
-// ----------------------------------------------------------------------
-void
-pylith::faults::FaultCohesiveTract::_initConstitutiveModel(void)
-{ // _initConstitutiveModel
- // :TODO: ADD STUFF HERE
-} // _initConstitutiveModel
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.hh 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.hh 2013-03-25 23:12:17 UTC (rev 21635)
@@ -51,11 +51,6 @@
virtual
void deallocate(void);
- /** Sets the spatial database for the inital tractions
- * @param dbs spatial database for initial tractions
- */
- void dbInitial(spatialdata::spatialdb::SpatialDB* dbs);
-
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
*
@@ -116,32 +111,6 @@
cellField(const char* name,
const topology::SolutionFields* fields =0);
- // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
-
- /** Calculate orientation at quadrature points.
- *
- * @param upDir Direction perpendicular to along-strike direction that is
- * not collinear with fault normal (usually "up" direction but could
- * be up-dip direction; applies to fault surfaces in 2-D and 3-D).
- */
- void _calcOrientation(const PylithScalar upDir[3]);
-
- /** Get initial tractions using a spatial database.
- */
- void _getInitialTractions(void);
-
- /** Setup fault constitutive model.
- */
- void _initConstitutiveModel(void);
-
- // PRIVATE MEMBERS ////////////////////////////////////////////////////
-private :
-
- /// Database for initial tractions.
- spatialdata::spatialdb::SpatialDB* _dbInitial;
-
-
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/LiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/LiuCosSlipFn.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/LiuCosSlipFn.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -23,6 +23,9 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/Stratum.hh" // USES Stratum
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -33,11 +36,6 @@
#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
-typedef pylith::topology::SubMesh::SieveMesh::label_sequence label_sequence;
-typedef pylith::topology::SubMesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::LiuCosSlipFn::LiuCosSlipFn(void) :
_slipTimeVertex(0),
@@ -70,42 +68,36 @@
// ----------------------------------------------------------------------
// Initialize slip time function.
void
-pylith::faults::LiuCosSlipFn::initialize(
- const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer,
- const PylithScalar originTime)
+pylith::faults::LiuCosSlipFn::initialize(const topology::SubMesh& faultMesh,
+ const spatialdata::units::Nondimensional& normalizer,
+ const PylithScalar originTime)
{ // initialize
- assert(0 != _dbFinalSlip);
- assert(0 != _dbSlipTime);
- assert(0 != _dbRiseTime);
+ assert(_dbFinalSlip);
+ assert(_dbSlipTime);
+ assert(_dbRiseTime);
- const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
- assert(0 != cs);
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();assert(cs);
const int spaceDim = cs->spaceDim();
const PylithScalar lengthScale = normalizer.lengthScale();
const PylithScalar timeScale = normalizer.timeScale();
// Get vertices in fault mesh
- DM dmMesh = faultMesh.dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = faultMesh.dmMesh();assert(dmMesh);
+ topology::Stratum depthStratum(dmMesh, topology::Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
- assert(dmMesh);
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-
delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
- assert(0 != _parameters);
+ assert(_parameters);
_parameters->add("final slip", "final_slip");
topology::Field<topology::SubMesh>& finalSlip = _parameters->get("final slip");
finalSlip.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
finalSlip.allocate();
finalSlip.scale(lengthScale);
finalSlip.vectorFieldType(topology::FieldBase::VECTOR);
- PetscSection finalSlipSection = finalSlip.petscSection();
- Vec finalSlipVec = finalSlip.localVector();
- PetscScalar *finalSlipArray;
- assert(finalSlipSection);assert(finalSlipVec);
+ topology::VecVisitorMesh finalSlipVisitor(finalSlip);
+ PetscScalar* finalSlipArray = finalSlipVisitor.localArray();
_parameters->add("slip time", "slip_time");
topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
@@ -113,10 +105,8 @@
slipTime.allocate();
slipTime.scale(timeScale);
slipTime.vectorFieldType(topology::FieldBase::SCALAR);
- PetscSection slipTimeSection = slipTime.petscSection();
- Vec slipTimeVec = slipTime.localVector();
- PetscScalar *slipTimeArray;
- assert(slipTimeSection);assert(slipTimeVec);
+ topology::VecVisitorMesh slipTimeVisitor(slipTime);
+ PetscScalar* slipTimeArray = slipTimeVisitor.localArray();
_parameters->add("rise time", "rise_time");
topology::Field<topology::SubMesh>& riseTime = _parameters->get("rise time");
@@ -124,10 +114,8 @@
riseTime.allocate();
riseTime.scale(timeScale);
riseTime.vectorFieldType(topology::FieldBase::SCALAR);
- PetscSection riseTimeSection = riseTime.petscSection();
- Vec riseTimeVec = riseTime.localVector();
- PetscScalar *riseTimeArray;
- assert(riseTimeSection);assert(riseTimeVec);
+ topology::VecVisitorMesh riseTimeVisitor(riseTime);
+ PetscScalar* riseTimeArray = riseTimeVisitor.localArray();
// Open databases and set query values
_dbFinalSlip->open();
@@ -165,36 +153,22 @@
// Get coordinates of vertices
scalar_array vCoordsGlobal(spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- PetscScalar *coordArray;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar* coordsArray = coordsVisitor.localArray();
_slipVertex.resize(spaceDim);
- err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(riseTimeVec, &riseTimeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt fsdof, fsoff, stdof, stoff, rtdof, rtoff, cdof, coff;
-
- err = PetscSectionGetDof(finalSlipSection, v, &fsdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(finalSlipSection, v, &fsoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(riseTimeSection, v, &rtdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(riseTimeSection, v, &rtoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(coordSection, v, &cdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(coordSection, v, &coff);CHECK_PETSC_ERROR(err);
- assert(cdof == spaceDim);
-
- for(PetscInt d = 0; d < cdof; ++d) {
- vCoordsGlobal[d] = coordArray[coff+d];
- }
+ // Dimensionalize coordinates
+ const PetscInt coff = coordsVisitor.sectionOffset(v);
+ assert(spaceDim == coordsVisitor.sectionDof(v));
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ vCoordsGlobal[d] = coordsArray[coff+d];
+ } // for
normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(), lengthScale);
+ // Final slip
+ const PetscInt fsoff = finalSlipVisitor.sectionOffset(v);
+ assert(spaceDim == finalSlipVisitor.sectionDof(v));
int err = _dbFinalSlip->query(&_slipVertex[0], _slipVertex.size(), &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -206,6 +180,9 @@
} // if
normalizer.nondimensionalize(&_slipVertex[0], _slipVertex.size(), lengthScale);
+ // Slip time
+ const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+ assert(1 == slipTimeVisitor.sectionDof(v));
err = _dbSlipTime->query(&_slipTimeVertex, 1, &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -219,6 +196,9 @@
// add origin time to rupture time
_slipTimeVertex += originTime;
+ // Rise time
+ const PetscInt rtoff = riseTimeVisitor.sectionOffset(v);
+ assert(1 == riseTimeVisitor.sectionDof(v));
err = _dbRiseTime->query(&_riseTimeVertex, 1, &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -230,16 +210,12 @@
} // if
normalizer.nondimensionalize(&_riseTimeVertex, 1, timeScale);
- for(PetscInt d = 0; d < fsdof; ++d) {
+ for(PetscInt d = 0; d < spaceDim; ++d) {
finalSlipArray[fsoff+d] = _slipVertex[d];
- }
+ } // for
slipTimeArray[stoff] = _slipTimeVertex;
riseTimeArray[stoff] = _riseTimeVertex;
} // for
- err = VecRestoreArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(riseTimeVec, &riseTimeArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
// Close databases
_dbFinalSlip->close();
@@ -251,59 +227,44 @@
// Get slip on fault surface at time t.
void
pylith::faults::LiuCosSlipFn::slip(topology::Field<topology::SubMesh>* slip,
- const PylithScalar t)
+ const PylithScalar t)
{ // slip
- assert(0 != slip);
- assert(0 != _parameters);
+ assert(slip);
+ assert(_parameters);
// Get vertices in fault mesh
- DM dmMesh = slip->mesh().dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = _parameters->mesh().dmMesh();assert(dmMesh);
+ topology::Stratum depthStratum(dmMesh, topology::Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
- assert(dmMesh);
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-
// Get sections
const topology::Field<topology::SubMesh>& finalSlip = _parameters->get("final slip");
- PetscSection finalSlipSection = finalSlip.petscSection();
- Vec finalSlipVec = finalSlip.localVector();
- PetscScalar *finalSlipArray;
- assert(finalSlipSection);assert(finalSlipVec);
+ topology::VecVisitorMesh finalSlipVisitor(finalSlip);
+ const PetscScalar* finalSlipArray = finalSlipVisitor.localArray();
+
const topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
- PetscSection slipTimeSection = slipTime.petscSection();
- Vec slipTimeVec = slipTime.localVector();
- PetscScalar *slipTimeArray;
- assert(slipTimeSection);assert(slipTimeVec);
+ topology::VecVisitorMesh slipTimeVisitor(slipTime);
+ const PetscScalar* slipTimeArray = slipTimeVisitor.localArray();
+
const topology::Field<topology::SubMesh>& riseTime = _parameters->get("rise time");
- PetscSection riseTimeSection = riseTime.petscSection();
- Vec riseTimeVec = riseTime.localVector();
- PetscScalar *riseTimeArray;
- assert(riseTimeSection);assert(riseTimeVec);
- PetscSection slipSection = slip->petscSection();
- Vec slipVec = slip->localVector();
- PetscScalar *slipArray;
- assert(slipSection);assert(slipVec);
+ topology::VecVisitorMesh riseTimeVisitor(riseTime);
+ const PetscScalar* riseTimeArray = riseTimeVisitor.localArray();
+ topology::VecVisitorMesh slipVisitor(*slip);
+ PetscScalar* slipArray = slipVisitor.localArray();
+
const int spaceDim = _slipVertex.size();
- err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(riseTimeVec, &riseTimeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt fsdof, fsoff, stdof, stoff, rtdof, rtoff, sdof, soff;
+ const PetscInt fsoff = finalSlipVisitor.sectionOffset(v);
+ const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+ const PetscInt rtoff = riseTimeVisitor.sectionOffset(v);
+ const PetscInt soff = slipVisitor.sectionOffset(v);
- err = PetscSectionGetDof(finalSlipSection, v, &fsdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(finalSlipSection, v, &fsoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(riseTimeSection, v, &rtdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(riseTimeSection, v, &rtoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipSection, v, &sdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipSection, v, &soff);CHECK_PETSC_ERROR(err);
- assert(fsdof == sdof);
- assert(stdof == 1);
- assert(rtdof == 1);
+ assert(spaceDim == finalSlipVisitor.sectionDof(v));
+ assert(1 == slipTimeVisitor.sectionDof(v));
+ assert(1 == riseTimeVisitor.sectionDof(v));
+ assert(spaceDim == slipVisitor.sectionDof(v));
PylithScalar finalSlipMag = 0.0;
for (int i=0; i < spaceDim; ++i)
@@ -314,9 +275,9 @@
const PylithScalar scale = finalSlipMag > 0.0 ? slip / finalSlipMag : 0.0;
// Update field
- for(PetscInt d = 0; d < sdof; ++d) {
+ for(PetscInt d = 0; d < spaceDim; ++d) {
slipArray[soff+d] += finalSlipArray[fsoff+d] * scale;
- }
+ } // for
} // for
PetscLogFlops((vEnd-vStart) * (2+28 + 3*_slipVertex.size()));
Deleted: short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -1,216 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard, U.S. Geological Survey
-// Charles A. Williams, GNS Science
-// Matthew G. Knepley, University of Chicago
-//
-// This code was developed as part of the Computational Infrastructure
-// for Geodynamics (http://geodynamics.org).
-//
-// Copyright (c) 2010-2012 University of California, Davis
-//
-// See COPYING for license information.
-//
-// ----------------------------------------------------------------------
-//
-
-#include <portinfo>
-
-#include "StaticPerturbation.hh" // implementation of object methods
-
-#include "pylith/topology/SubMesh.hh" // USES SubMesh
-#include "pylith/topology/Fields.hh" // USES Fields
-#include "pylith/topology/Field.hh" // USES Field
-
-#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
-#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
-#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-
-#include <cassert> // USES assert()
-#include <sstream> // USES std::ostringstream
-#include <stdexcept> // USES std::runtime_error
-
-// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
-typedef pylith::topology::SubMesh::SieveMesh::label_sequence label_sequence;
-typedef pylith::topology::SubMesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
-// Default constructor.
-pylith::faults::StaticPerturbation::StaticPerturbation(void) :
- _dbAmplitude(0)
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor.
-pylith::faults::StaticPerturbation::~StaticPerturbation(void)
-{ // destructor
- deallocate();
-} // destructor
-
-// ----------------------------------------------------------------------
-// Deallocate PETSc and local data structures.
-void
-pylith::faults::StaticPerturbation::deallocate(void)
-{ // deallocate
- TractPerturbation::deallocate();
-
- _dbAmplitude = 0; // :TODO: Use shared pointer
-} // deallocate
-
-// ----------------------------------------------------------------------
-// Set spatial database for amplitude of traction.
-void
-pylith::faults::StaticPerturbation::dbAmplitude(spatialdata::spatialdb::SpatialDB* const db) {
- _dbAmplitude = db;
-} // dbAmplitude
-
-// ----------------------------------------------------------------------
-// Initialize traction perturbation function.
-void
-pylith::faults::StaticPerturbation::initialize(
- const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer)
-{ // initialize
- assert(_dbAmplitude);
-
- const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
- assert(cs);
- const int spaceDim = cs->spaceDim();
-
- const PylithScalar lengthScale = normalizer.lengthScale();
- const PylithScalar timeScale = normalizer.timeScale();
- const PylithScalar pressureScale = normalizer.pressureScale();
-
- // Get vertices in fault mesh
- const ALE::Obj<SieveMesh>& sieveMesh = faultMesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<label_sequence>& vertices = sieveMesh->depthStratum(0);
- assert(!vertices.isNull());
- const label_sequence::iterator verticesBegin = vertices->begin();
- const label_sequence::iterator verticesEnd = vertices->end();
-
- delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
- assert(_parameters);
- _parameters->add("amplitude", "amplitude");
-
- topology::Field<topology::SubMesh>& amplitude = _parameters->get("amplitude");
- amplitude.newSection(vertices, spaceDim);
- amplitude.allocate();
- amplitude.scale(pressureScale);
- amplitude.vectorFieldType(topology::FieldBase::VECTOR);
- const ALE::Obj<RealSection>& amplitudeSection = amplitude.section();
- assert(!amplitudeSection.isNull());
-
- // Open databases and set query values
- _dbAmplitude->open();
- switch (spaceDim)
- { // switch
- case 1 : {
- const char* tractionValues[1] = {"traction-normal"};
- _dbAmplitude->queryVals(tractionValues, 1);
- break;
- } // case 1
- case 2 : {
- const char* tractionValues[2] = {"traction-shear", "traction-normal"};
- _dbAmplitude->queryVals(tractionValues, 2);
- break;
- } // case 2
- case 3 : {
- const char* tractionValues[3] = {"traction-shear-leftlateral",
- "traction-shear-updip",
- "traction-normal"};
- _dbAmplitude->queryVals(tractionValues, 3);
- break;
- } // case 3
- default :
- std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
- assert(0);
- throw std::logic_error("Bad spatial dimension in StaticPerturbation.");
- } // switch
-
- // Get coordinates of vertices
- const ALE::Obj<RealSection>& coordinates = sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
-
- scalar_array tractionVertex(spaceDim);
- scalar_array vCoordsGlobal(spaceDim);
- for (label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- coordinates->restrictPoint(*v_iter, &vCoordsGlobal[0], vCoordsGlobal.size());
- normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(), lengthScale);
-
- int err = _dbAmplitude->query(&tractionVertex[0], tractionVertex.size(),
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
- if (err) {
- std::ostringstream msg;
- msg << "Could not find traction at (";
- for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoordsGlobal[i];
- msg << ") using spatial database " << _dbAmplitude->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- normalizer.nondimensionalize(&tractionVertex[0], tractionVertex.size(), pressureScale);
-
- assert(spaceDim == amplitudeSection->getFiberDimension(*v_iter));
- amplitudeSection->updatePoint(*v_iter, &tractionVertex[0]);
- } // for
-
- // Close databases
- _dbAmplitude->close();
-} // initialize
-
-// ----------------------------------------------------------------------
-// Get traction on fault surface at time t.
-void
-pylith::faults::StaticPerturbation::traction(topology::Field<topology::SubMesh>* tractionField,
- const PylithScalar t)
-{ // traction
- assert(tractionField);
- assert(_parameters);
-
- const spatialdata::geocoords::CoordSys* cs = tractionField->mesh().coordsys();
- assert(cs);
- const int spaceDim = cs->spaceDim();
-
- // Get vertices in fault mesh
- const ALE::Obj<SieveMesh>& sieveMesh = tractionField->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<label_sequence>& vertices = sieveMesh->depthStratum(0);
- assert(!vertices.isNull());
- const label_sequence::iterator verticesBegin = vertices->begin();
- const label_sequence::iterator verticesEnd = vertices->end();
-
- // Get sections
- const ALE::Obj<RealSection>& amplitudeSection = _parameters->get("amplitude").section();
- assert(!amplitudeSection.isNull());
- const ALE::Obj<RealSection>& tractionSection = tractionField->section();
- assert(!tractionSection.isNull());
-
- for (label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- assert(spaceDim == amplitudeSection->getFiberDimension(*v_iter));
- const PylithScalar* amplitudeVertex = amplitudeSection->restrictPoint(*v_iter);
-
- // Update field
- assert(spaceDim == tractionSection->getFiberDimension(*v_iter));
- tractionSection->updateAddPoint(*v_iter, &litudeVertex[0]);
- } // for
-
-} // traction
-
-// ----------------------------------------------------------------------
-// Get traction amplitude..
-const pylith::topology::Field<pylith::topology::SubMesh>&
-pylith::faults::StaticPerturbation::amplitude(void)
-{ // amplitude
- return _parameters->get("amplitude");
-} // amplitude
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.hh 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.hh 2013-03-25 23:12:17 UTC (rev 21635)
@@ -1,102 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard, U.S. Geological Survey
-// Charles A. Williams, GNS Science
-// Matthew G. Knepley, University of Chicago
-//
-// This code was developed as part of the Computational Infrastructure
-// for Geodynamics (http://geodynamics.org).
-//
-// Copyright (c) 2010-2012 University of California, Davis
-//
-// See COPYING for license information.
-//
-// ----------------------------------------------------------------------
-//
-
-/** @file libsrc/faults/StaticPerturbation.hh
- *
- * @brief C++ implementation of a static perturbation in tractions.
- */
-
-#if !defined(pylith_faults_staticperturbation_hh)
-#define pylith_faults_staticperturbation_hh
-
-// Include directives ---------------------------------------------------
-#include "TractPerturbation.hh"
-
-#include "pylith/utils/array.hh" // HASA scalar_array
-
-#include "spatialdata/spatialdb/spatialdbfwd.hh"
-
-// StaticPerturbation -----------------------------------------------------------
-/**
- * @brief C++ implementation of a static perturbation in tractions.
- *
- * T = F(x)
-*/
-class pylith::faults::StaticPerturbation : public TractPerturbation
-{ // class StaticPerturbation
- friend class TestStaticPerturbation; // unit testing
-
-// PUBLIC METHODS ///////////////////////////////////////////////////////
-public :
-
- /// Default constructor.
- StaticPerturbation(void);
-
- /// Destructor.
- ~StaticPerturbation(void);
-
- /// Deallocate PETSc and local data structures.
- virtual
- void deallocate(void);
-
- /** Set spatial database for traction amplitude.
- *
- * @param db Spatial database
- */
- void dbAmplitude(spatialdata::spatialdb::SpatialDB* const db);
-
- /** Initialize static perturbation.
- *
- * @param faultMesh Finite-element mesh of fault.
- * @param normalizer Nondimensionalization of scales.
- */
- void initialize(const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer);
-
- /** Get traction on fault surface at time t.
- *
- * @param tractionField Traction field over fault surface.
- * @param t Time t.
- */
- void traction(topology::Field<topology::SubMesh>* const tractionField,
- const PylithScalar t);
-
- /** Get amplitude of traction perturbation.
- *
- * @returns Final slip.
- */
- const topology::Field<topology::SubMesh>& amplitude(void);
-
-// NOT IMPLEMENTED //////////////////////////////////////////////////////
-private :
-
- StaticPerturbation(const StaticPerturbation&); ///< Not implemented.
- const StaticPerturbation& operator=(const StaticPerturbation&); ///< Not implemented
-
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
-
- /// Spatial database for traction perturbation.
- spatialdata::spatialdb::SpatialDB* _dbAmplitude;
-
-}; // class StaticPerturbation
-
-#endif // pylith_faults_staticperturbation_hh
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.icc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/StaticPerturbation.icc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -1,38 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard, U.S. Geological Survey
-// Charles A. Williams, GNS Science
-// Matthew G. Knepley, University of Chicago
-//
-// This code was developed as part of the Computational Infrastructure
-// for Geodynamics (http://geodynamics.org).
-//
-// Copyright (c) 2010-2011 University of California, Davis
-//
-// See COPYING for license information.
-//
-// ----------------------------------------------------------------------
-//
-
-#if !defined(pylith_faults_stepslipfn_hh)
-#error "StepSlipFn.icc can only be included from StepSlipFn.hh"
-#endif
-
-// Set spatial database for final slip.
-inline
-void
-pylith::faults::StepSlipFn::dbFinalSlip(spatialdata::spatialdb::SpatialDB* const db) {
- _dbFinalSlip = db;
-} // dbFinalSlip
-
-// Set spatial database for slip initiation time.
-inline
-void
-pylith::faults::StepSlipFn::dbSlipTime(spatialdata::spatialdb::SpatialDB* const db) {
- _dbSlipTime = db;
-} // dbSlipTime
-
-
-// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/StepSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/StepSlipFn.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/StepSlipFn.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -23,6 +23,9 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/Stratum.hh" // USES Stratum
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -33,11 +36,6 @@
#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
-typedef pylith::topology::SubMesh::SieveMesh::label_sequence label_sequence;
-typedef pylith::topology::SubMesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::StepSlipFn::StepSlipFn(void) :
_slipTimeVertex(0),
@@ -67,31 +65,27 @@
// ----------------------------------------------------------------------
// Initialize slip time function.
void
-pylith::faults::StepSlipFn::initialize(
- const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer,
- const PylithScalar originTime)
+pylith::faults::StepSlipFn::initialize(const topology::SubMesh& faultMesh,
+ const spatialdata::units::Nondimensional& normalizer,
+ const PylithScalar originTime)
{ // initialize
- assert(0 != _dbFinalSlip);
- assert(0 != _dbSlipTime);
+ assert(_dbFinalSlip);
+ assert(_dbSlipTime);
- const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
- assert(0 != cs);
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();assert(cs);
const int spaceDim = cs->spaceDim();
const PylithScalar lengthScale = normalizer.lengthScale();
const PylithScalar timeScale = normalizer.timeScale();
// Get vertices in fault mesh
- DM dmMesh = faultMesh.dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = faultMesh.dmMesh();assert(dmMesh);
+ topology::Stratum depthStratum(dmMesh, topology::Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
- assert(dmMesh);
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-
delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
- assert(0 != _parameters);
+ assert(_parameters);
_parameters->add("final slip", "final_slip");
topology::Field<topology::SubMesh>& finalSlip = _parameters->get("final slip");
@@ -99,10 +93,8 @@
finalSlip.allocate();
finalSlip.scale(lengthScale);
finalSlip.vectorFieldType(topology::FieldBase::VECTOR);
- PetscSection finalSlipSection = finalSlip.petscSection();
- Vec finalSlipVec = finalSlip.localVector();
- PetscScalar *finalSlipArray;
- assert(finalSlipSection);assert(finalSlipVec);
+ topology::VecVisitorMesh finalSlipVisitor(finalSlip);
+ PetscScalar* finalSlipArray = finalSlipVisitor.localArray();
_parameters->add("slip time", "slip_time");
topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
@@ -110,10 +102,8 @@
slipTime.allocate();
slipTime.scale(timeScale);
slipTime.vectorFieldType(topology::FieldBase::SCALAR);
- PetscSection slipTimeSection = slipTime.petscSection();
- Vec slipTimeVec = slipTime.localVector();
- PetscScalar *slipTimeArray;
- assert(slipTimeSection);assert(slipTimeVec);
+ topology::VecVisitorMesh slipTimeVisitor(slipTime);
+ PetscScalar* slipTimeArray = slipTimeVisitor.localArray();
// Open databases and set query values
_dbFinalSlip->open();
@@ -147,33 +137,22 @@
// Get coordinates of vertices
scalar_array vCoordsGlobal(spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- PetscScalar *coordArray;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar* coordsArray = coordsVisitor.localArray();
_slipVertex.resize(spaceDim);
- err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt fsdof, fsoff, stdof, stoff, cdof, coff;
-
- err = PetscSectionGetDof(finalSlipSection, v, &fsdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(finalSlipSection, v, &fsoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(coordSection, v, &cdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(coordSection, v, &coff);CHECK_PETSC_ERROR(err);
- assert(cdof == spaceDim);
-
- for(PetscInt d = 0; d < cdof; ++d) {
- vCoordsGlobal[d] = coordArray[coff+d];
- }
+ // Dimensionalize coordinates
+ const PetscInt coff = coordsVisitor.sectionOffset(v);
+ assert(spaceDim == coordsVisitor.sectionDof(v));
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ vCoordsGlobal[d] = coordsArray[coff+d];
+ } // for
normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(), lengthScale);
-
+
+ // Final slip
+ const PetscInt fsoff = finalSlipVisitor.sectionOffset(v);
+ assert(spaceDim == finalSlipVisitor.sectionDof(v));
int err = _dbFinalSlip->query(&_slipVertex[0], _slipVertex.size(), &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -185,6 +164,9 @@
} // if
normalizer.nondimensionalize(&_slipVertex[0], _slipVertex.size(), lengthScale);
+ // Slip time
+ const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+ assert(1 == slipTimeVisitor.sectionDof(v));
err = _dbSlipTime->query(&_slipTimeVertex, 1, &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -198,14 +180,11 @@
// add origin time to rupture time
_slipTimeVertex += originTime;
- for(PetscInt d = 0; d < fsdof; ++d) {
+ for(PetscInt d = 0; d < spaceDim; ++d) {
finalSlipArray[fsoff+d] = _slipVertex[d];
- }
+ } // for
slipTimeArray[stoff] = _slipTimeVertex;
} // for
- err = VecRestoreArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
// Close databases
_dbFinalSlip->close();
@@ -218,58 +197,44 @@
pylith::faults::StepSlipFn::slip(topology::Field<topology::SubMesh>* slip,
const PylithScalar t)
{ // slip
- assert(0 != slip);
- assert(0 != _parameters);
+ assert(slip);
+ assert(_parameters);
// Get vertices in fault mesh
- DM dmMesh = slip->mesh().dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = _parameters->mesh().dmMesh();assert(dmMesh);
+ topology::Stratum depthStratum(dmMesh, topology::Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
- assert(dmMesh);
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-
// Get sections
const topology::Field<topology::SubMesh>& finalSlip = _parameters->get("final slip");
- const topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
- PetscSection finalSlipSection = finalSlip.petscSection();
- Vec finalSlipVec = finalSlip.localVector();
- PetscScalar *finalSlipArray;
- assert(finalSlipSection);assert(finalSlipVec);
- PetscSection slipTimeSection = slipTime.petscSection();
- Vec slipTimeVec = slipTime.localVector();
- PetscScalar *slipTimeArray;
- assert(slipTimeSection);assert(slipTimeVec);
- PetscSection slipSection = slip->petscSection();
- Vec slipVec = slip->localVector();
- PetscScalar *slipArray;
- assert(slipSection);assert(slipVec);
+ topology::VecVisitorMesh finalSlipVisitor(finalSlip);
+ const PetscScalar* finalSlipArray = finalSlipVisitor.localArray();
- err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+ const topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
+ topology::VecVisitorMesh slipTimeVisitor(slipTime);
+ const PetscScalar* slipTimeArray = slipTimeVisitor.localArray();
+
+ topology::VecVisitorMesh slipVisitor(*slip);
+ PetscScalar* slipArray = slipVisitor.localArray();
+
+ const int spaceDim = _slipVertex.size();
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt fsdof, fsoff, stdof, stoff, sdof, soff;
+ const PetscInt fsoff = finalSlipVisitor.sectionOffset(v);
+ const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+ const PetscInt soff = slipVisitor.sectionOffset(v);
- err = PetscSectionGetDof(finalSlipSection, v, &fsdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(finalSlipSection, v, &fsoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipSection, v, &sdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipSection, v, &soff);CHECK_PETSC_ERROR(err);
- assert(stdof == 1);
- assert(fsdof == sdof);
+ assert(spaceDim == finalSlipVisitor.sectionDof(v));
+ assert(1 == slipTimeVisitor.sectionDof(v));
+ assert(spaceDim == slipVisitor.sectionDof(v));
const PylithScalar relTime = t - slipTimeArray[stoff];
if (relTime >= 0.0) {
- for(PetscInt d = 0; d < sdof; ++d) {
+ for(PetscInt d = 0; d < spaceDim; ++d) {
slipArray[soff+d] += finalSlipArray[fsoff+d];
- }
- }
+ } // for
+ } // if
} // for
- err = VecRestoreArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
PetscLogFlops((vEnd-vStart) * 1);
} // slip
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/TimeHistorySlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/TimeHistorySlipFn.cc 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/TimeHistorySlipFn.cc 2013-03-25 23:12:17 UTC (rev 21635)
@@ -23,6 +23,9 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/Stratum.hh" // USES Stratum
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/spatialdb/TimeHistory.hh" // USES TimeHistory
@@ -34,11 +37,6 @@
#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
-typedef pylith::topology::SubMesh::SieveMesh::label_sequence label_sequence;
-typedef pylith::topology::SubMesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::TimeHistorySlipFn::TimeHistorySlipFn(void) :
_slipTimeVertex(0),
@@ -65,7 +63,7 @@
_dbAmplitude = 0; // :TODO: Use shared pointer
_dbSlipTime = 0; // :TODO: Use shared pointer
- if (0 != _dbTimeHistory)
+ if (_dbTimeHistory)
_dbTimeHistory->close();
_dbTimeHistory = 0; // :TODO: Use shared pointer
} // deallocate
@@ -73,31 +71,27 @@
// ----------------------------------------------------------------------
// Initialize slip time function.
void
-pylith::faults::TimeHistorySlipFn::initialize(
- const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer,
- const PylithScalar originTime)
+pylith::faults::TimeHistorySlipFn::initialize(const topology::SubMesh& faultMesh,
+ const spatialdata::units::Nondimensional& normalizer,
+ const PylithScalar originTime)
{ // initialize
- assert(0 != _dbAmplitude);
- assert(0 != _dbSlipTime);
+ assert(_dbAmplitude);
+ assert(_dbSlipTime);
- const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
- assert(0 != cs);
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();assert(cs);
const int spaceDim = cs->spaceDim();
const PylithScalar lengthScale = normalizer.lengthScale();
const PylithScalar timeScale = normalizer.timeScale();
// Get vertices in fault mesh
- DM dmMesh = faultMesh.dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = faultMesh.dmMesh();assert(dmMesh);
+ topology::Stratum depthStratum(dmMesh, topology::Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
- assert(dmMesh);
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-
delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
- assert(0 != _parameters);
+ assert(_parameters);
_parameters->add("slip amplitude", "slip_amplitude");
topology::Field<topology::SubMesh>& slipAmplitude = _parameters->get("slip amplitude");
@@ -105,10 +99,8 @@
slipAmplitude.allocate();
slipAmplitude.scale(lengthScale);
slipAmplitude.vectorFieldType(topology::FieldBase::VECTOR);
- PetscSection finalSlipSection = slipAmplitude.petscSection();
- Vec finalSlipVec = slipAmplitude.localVector();
- PetscScalar *finalSlipArray;
- assert(finalSlipSection);assert(finalSlipVec);
+ topology::VecVisitorMesh slipAmplitudeVisitor(slipAmplitude);
+ PetscScalar* slipAmplitudeArray = slipAmplitudeVisitor.localArray();
_parameters->add("slip time", "slip_time");
topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
@@ -116,10 +108,8 @@
slipTime.allocate();
slipTime.scale(timeScale);
slipTime.vectorFieldType(topology::FieldBase::SCALAR);
- PetscSection slipTimeSection = slipTime.petscSection();
- Vec slipTimeVec = slipTime.localVector();
- PetscScalar *slipTimeArray;
- assert(slipTimeSection);assert(slipTimeVec);
+ topology::VecVisitorMesh slipTimeVisitor(slipTime);
+ PetscScalar* slipTimeArray = slipTimeVisitor.localArray();
// Open databases and set query values
_dbAmplitude->open();
@@ -153,33 +143,22 @@
// Get coordinates of vertices
scalar_array vCoordsGlobal(spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- PetscScalar *coordArray;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar* coordsArray = coordsVisitor.localArray();
_slipVertex.resize(spaceDim);
- err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt fsdof, fsoff, stdof, stoff, cdof, coff;
-
- err = PetscSectionGetDof(finalSlipSection, v, &fsdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(finalSlipSection, v, &fsoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(coordSection, v, &cdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(coordSection, v, &coff);CHECK_PETSC_ERROR(err);
- assert(cdof == spaceDim);
-
- for(PetscInt d = 0; d < cdof; ++d) {
- vCoordsGlobal[d] = coordArray[coff+d];
- }
+ // Dimensionalize coordinates
+ const PetscInt coff = coordsVisitor.sectionOffset(v);
+ assert(spaceDim == coordsVisitor.sectionDof(v));
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ vCoordsGlobal[d] = coordsArray[coff+d];
+ } // for
normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(), lengthScale);
-
+
+ // Slip amplitude
+ const PetscInt saoff = slipAmplitudeVisitor.sectionOffset(v);
+ assert(spaceDim == slipAmplitudeVisitor.sectionDof(v));
int err = _dbAmplitude->query(&_slipVertex[0], _slipVertex.size(), &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -191,8 +170,10 @@
} // if
normalizer.nondimensionalize(&_slipVertex[0], _slipVertex.size(), lengthScale);
- err = _dbSlipTime->query(&_slipTimeVertex, 1,
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+ // Slip time
+ const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+ assert(1 == slipTimeVisitor.sectionDof(v));
+ err = _dbSlipTime->query(&_slipTimeVertex, 1, &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
msg << "Could not find slip initiation time at (";
@@ -205,9 +186,9 @@
// add origin time to rupture time
_slipTimeVertex += originTime;
- for(PetscInt d = 0; d < fsdof; ++d) {
- finalSlipArray[fsoff+d] = _slipVertex[d];
- }
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ slipAmplitudeArray[saoff+d] = _slipVertex[d];
+ } // for
slipTimeArray[stoff] = _slipTimeVertex;
} // for
@@ -224,57 +205,46 @@
// Get slip on fault surface at time t.
void
pylith::faults::TimeHistorySlipFn::slip(topology::Field<topology::SubMesh>* slip,
- const PylithScalar t)
+ const PylithScalar t)
{ // slip
- assert(0 != slip);
- assert(0 != _parameters);
- assert(0 != _dbTimeHistory);
+ assert(slip);
+ assert(_parameters);
+ assert(_dbTimeHistory);
// Get vertices in fault mesh
- DM dmMesh = slip->mesh().dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = _parameters->mesh().dmMesh();assert(dmMesh);
+ topology::Stratum depthStratum(dmMesh, topology::Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
- assert(dmMesh);
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
-
// Get sections
- const topology::Field<topology::SubMesh>& finalSlip = _parameters->get("slip amplitude");
- PetscSection finalSlipSection = finalSlip.petscSection();
- Vec finalSlipVec = finalSlip.localVector();
- PetscScalar *finalSlipArray;
- assert(finalSlipSection);assert(finalSlipVec);
+ const topology::Field<topology::SubMesh>& slipAmplitude = _parameters->get("slip amplitude");
+ topology::VecVisitorMesh slipAmplitudeVisitor(slipAmplitude);
+ const PetscScalar* slipAmplitudeArray = slipAmplitudeVisitor.localArray();
+
const topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
- PetscSection slipTimeSection = slipTime.petscSection();
- Vec slipTimeVec = slipTime.localVector();
- PetscScalar *slipTimeArray;
- assert(slipTimeSection);assert(slipTimeVec);
- PetscSection slipSection = slip->petscSection();
- Vec slipVec = slip->localVector();
- PetscScalar *slipArray;
- assert(slipSection);assert(slipVec);
+ topology::VecVisitorMesh slipTimeVisitor(slipTime);
+ const PetscScalar* slipTimeArray = slipTimeVisitor.localArray();
+ topology::VecVisitorMesh slipVisitor(*slip);
+ PetscScalar* slipArray = slipVisitor.localArray();
+
+ const int spaceDim = _slipVertex.size();
PylithScalar amplitude = 0.0;
- err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt fsdof, fsoff, stdof, stoff, sdof, soff;
+ const PetscInt saoff = slipAmplitudeVisitor.sectionOffset(v);
+ const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+ const PetscInt soff = slipVisitor.sectionOffset(v);
- err = PetscSectionGetDof(finalSlipSection, v, &fsdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(finalSlipSection, v, &fsoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipTimeSection, v, &stdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipTimeSection, v, &stoff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(slipSection, v, &sdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(slipSection, v, &soff);CHECK_PETSC_ERROR(err);
- assert(fsdof == sdof);
- assert(stdof == 1);
+ assert(spaceDim == slipAmplitudeVisitor.sectionDof(v));
+ assert(1 == slipTimeVisitor.sectionDof(v));
+ assert(spaceDim == slipVisitor.sectionDof(v));
PylithScalar relTime = t - slipTimeArray[stoff];
if (relTime >= 0.0) {
relTime *= _timeScale;
const int err = _dbTimeHistory->query(&litude, relTime);
- if (0 != err) {
+ if (err) {
std::ostringstream msg;
msg << "Error querying for time '" << relTime
<< "' in time history database "
@@ -282,9 +252,9 @@
throw std::runtime_error(msg.str());
} // if
- for(PetscInt d = 0; d < sdof; ++d) {
- slipArray[soff+d] += finalSlipArray[fsoff+d] * amplitude;
- }
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ slipArray[soff+d] += slipAmplitudeArray[saoff+d] * amplitude;
+ } // for
} // else
} // for
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/faultsfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/faultsfwd.hh 2013-03-25 22:59:33 UTC (rev 21634)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/faultsfwd.hh 2013-03-25 23:12:17 UTC (rev 21635)
@@ -50,8 +50,6 @@
class Nucleator;
class TractPerturbation;
- class StaticPerturbation;
- class SpaceTimePerturbation;
class TopologyOps;
template<typename Sieve, typename Renumbering> class ReplaceVisitor;
More information about the CIG-COMMITS
mailing list