[cig-commits] r21608 - in short/3D/PyLith/trunk/libsrc/pylith: bc faults
brad at geodynamics.org
brad at geodynamics.org
Thu Mar 21 18:48:52 PDT 2013
Author: brad
Date: 2013-03-21 18:48:52 -0700 (Thu, 21 Mar 2013)
New Revision: 21608
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc
Log:
Code cleanup. Update to use visitors.
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2013-03-22 00:32:36 UTC (rev 21607)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2013-03-22 01:48:52 UTC (rev 21608)
@@ -26,6 +26,7 @@
#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
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -227,12 +228,11 @@
assert(_parameters);
const topology::Mesh& mesh = _parameters->mesh();
- const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
- assert(cs);
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();assert(cs);
const int spaceDim = cs->spaceDim();
- const PylithScalar lengthScale = _getNormalizer().lengthScale();
- PetscErrorCode err = 0;
+ const spatialdata::units::Nondimensional& normalizer = _getNormalizer();
+ const PylithScalar lengthScale = normalizer.lengthScale();
scalar_array coordsVertex(spaceDim);
PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
@@ -252,7 +252,8 @@
for (PetscInt d = 0; d < spaceDim; ++d) {
coordsVertex[d] = coordArray[coff+d];
} // for
- _getNormalizer().dimensionalize(&coordsVertex[0], coordsVertex.size(), lengthScale);
+ normalizer.dimensionalize(&coordsVertex[0], coordsVertex.size(), lengthScale);
+
int err = db->query(&valueVertex[0], valueVertex.size(), &coordsVertex[0], coordsVertex.size(), cs);
if (err) {
std::ostringstream msg;
@@ -262,12 +263,12 @@
msg << ") using spatial database " << db->label() << ".";
throw std::runtime_error(msg.str());
} // if
- _getNormalizer().nondimensionalize(&valueVertex[0], valueVertex.size(), scale);
+ normalizer.nondimensionalize(&valueVertex[0], valueVertex.size(), scale);
// Update section
const PetscInt off = parametersVisitor.sectionOffset(_points[iPoint]);
- const PetscInt dof = parametersVisitor.sectionDof(_points[iPoint]);
- for(int i = 0; i < dof; ++i) {
+ assert(querySize == parametersVisitor.sectionDof(_points[iPoint]));
+ for(int i = 0; i < querySize; ++i) {
parametersArray[off+i] = valueVertex[i];
} // for
} // for
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc 2013-03-22 00:32:36 UTC (rev 21607)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc 2013-03-22 01:48:52 UTC (rev 21608)
@@ -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
@@ -65,17 +68,15 @@
// ----------------------------------------------------------------------
// Initialize slip time function.
void
-pylith::faults::BruneSlipFn::initialize(
- const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer,
- const PylithScalar originTime)
+pylith::faults::BruneSlipFn::initialize(const topology::SubMesh& faultMesh,
+ const spatialdata::units::Nondimensional& normalizer,
+ const PylithScalar originTime)
{ // initialize
assert(_dbFinalSlip);
assert(_dbSlipTime);
assert(_dbRiseTime);
- 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();
@@ -83,9 +84,9 @@
// 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);
@@ -95,9 +96,8 @@
finalSlip.allocate();
finalSlip.scale(lengthScale);
finalSlip.vectorFieldType(topology::FieldBase::VECTOR);
- PetscSection finalSlipSection = finalSlip.petscSection();assert(finalSlipSection);
- PetscVec finalSlipVec = finalSlip.localVector();assert(finalSlipVec);
- PetscScalar *finalSlipArray = NULL;
+ topology::VecVisitorMesh finalSlipVisitor(finalSlip);
+ PetscScalar* finalSlipArray = finalSlipVisitor.localArray();
_parameters->add("slip time", "slip_time");
topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
@@ -105,9 +105,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();
_parameters->add("rise time", "rise_time");
topology::Field<topology::SubMesh>& riseTime = _parameters->get("rise time");
@@ -115,9 +114,8 @@
riseTime.allocate();
riseTime.scale(timeScale);
riseTime.vectorFieldType(topology::FieldBase::SCALAR);
- PetscSection riseTimeSection = riseTime.petscSection();assert(riseTimeSection);
- PetscVec riseTimeVec = riseTime.localVector();assert(riseTimeVec);
- PetscScalar *riseTimeArray = NULL;
+ topology::VecVisitorMesh riseTimeVisitor(riseTime);
+ PetscScalar* riseTimeArray = riseTimeVisitor.localArray();
// Open databases and set query values
_dbFinalSlip->open();
@@ -141,7 +139,7 @@
} // case 3
default :
std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
- assert(0);
+ assert(false);
throw std::logic_error("Bad spatial dimension in BruneSlipFn.");
} // switch
@@ -155,35 +153,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);
- err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
+ 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);
-
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;
@@ -195,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;
@@ -208,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;
@@ -219,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];
}
slipTimeArray[stoff] = _slipTimeVertex;
- riseTimeArray[stoff] = _riseTimeVertex;
+ riseTimeArray[rtoff] = _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();
@@ -246,50 +233,40 @@
assert(_parameters);
// Get vertices in fault mesh
- PetscDM dmMesh = slip->mesh().dmMesh();assert(dmMesh);
- PetscInt vStart, vEnd;
+ 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();
+
PetscErrorCode err;
- 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();assert(finalSlipSection);
- PetscVec finalSlipVec = finalSlip.localVector();assert(finalSlipVec);
- PetscScalar *finalSlipArray = NULL;
- err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
+ topology::VecVisitorMesh finalSlipVisitor(finalSlip);
+ const PetscScalar* finalSlipArray = finalSlipVisitor.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();
const topology::Field<topology::SubMesh>& riseTime = _parameters->get("rise time");
- PetscSection riseTimeSection = riseTime.petscSection();assert(riseTimeSection);
- PetscVec riseTimeVec = riseTime.localVector();assert(riseTimeVec);
- PetscScalar *riseTimeArray = NULL;
- err = VecGetArray(riseTimeVec, &riseTimeArray);CHECK_PETSC_ERROR(err);
+ topology::VecVisitorMesh riseTimeVisitor(riseTime);
+ const PetscScalar* riseTimeArray = riseTimeVisitor.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 = _slipVertex.size();
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt fsdof, fsoff, stdof, stoff, rtdof, rtoff, sdof, soff;
- 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);
+ const PetscInt fsoff = finalSlipVisitor.sectionOffset(v);
+ const PetscInt stoff = slipTimeVisitor.sectionOffset(v);
+ const PetscInt rtoff = riseTimeVisitor.sectionOffset(v);
+ const PetscInt soff = slipVisitor.sectionOffset(v);
- 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)
@@ -300,14 +277,10 @@
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
- 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(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
PetscLogFlops((vEnd-vStart) * (2+8 + 3*spaceDim));
} // slip
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc 2013-03-22 00:32:36 UTC (rev 21607)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc 2013-03-22 01:48:52 UTC (rev 21608)
@@ -25,6 +25,9 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Fields.hh" // HOLDSA 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
@@ -252,130 +255,105 @@
const PylithScalar timeScale = _timeScale;
// Get vertices.
- DM dmMesh = _parameters->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);
-
- const spatialdata::geocoords::CoordSys* cs = _parameters->mesh().coordsys();
- assert(cs);
+ const spatialdata::geocoords::CoordSys* cs = _parameters->mesh().coordsys();assert(cs);
const int spaceDim = cs->spaceDim();
- PetscSection initialSection, rateSection, rateTimeSection, changeSection, changeTimeSection;
- Vec initialVec, rateVec, rateTimeVec, changeVec, changeTimeVec;
- PetscScalar *valuesArray, *initialArray, *rateArray, *rateTimeArray, *changeArray, *changeTimeArray;
+ topology::Field<topology::SubMesh>& valueField = _parameters->get("value");
+ topology::VecVisitorMesh valueVisitor(valueField);
+ PetscScalar* valueArray = valueVisitor.localArray();
+
+ topology::Field<topology::SubMesh>* initialField = (_dbInitial) ? &_parameters->get("initial") : 0;
+ topology::VecVisitorMesh* initialVisitor = (initialField) ? new topology::VecVisitorMesh(*initialField) : 0;
+ PetscScalar* initialArray = (initialVisitor) ? initialVisitor->localArray() : NULL;
- PetscSection valuesSection = _parameters->get("value").petscSection();
- Vec valuesVec = _parameters->get("value").localVector();
- assert(valuesSection);assert(valuesVec);
- err = VecGetArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
- if (_dbInitial) {
- initialSection = _parameters->get("initial").petscSection();
- initialVec = _parameters->get("initial").localVector();
- assert(initialSection);assert(initialVec);
- err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
- }
- if (_dbRate) {
- rateSection = _parameters->get("rate").petscSection();
- rateTimeSection = _parameters->get("rate time").petscSection();
- rateVec = _parameters->get("rate").localVector();
- rateTimeVec = _parameters->get("rate time").localVector();
- assert(rateSection);assert(rateTimeSection);assert(rateVec);assert(rateTimeVec);
- err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
- }
- if (_dbChange) {
- changeSection = _parameters->get("change").petscSection();
- changeTimeSection = _parameters->get("change time").petscSection();
- changeVec = _parameters->get("change").localVector();
- changeTimeVec = _parameters->get("change time").localVector();
- assert(changeSection);assert(changeTimeSection);assert(changeVec);assert(changeTimeVec);
- err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
- }
+ topology::Field<topology::SubMesh>* rateField = (_dbRate) ? &_parameters->get("rate") : 0;
+ topology::VecVisitorMesh* rateVisitor = (rateField) ? new topology::VecVisitorMesh(*rateField) : 0;
+ PetscScalar* rateArray = (rateVisitor) ? rateVisitor->localArray() : NULL;
- for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt vdof, voff;
+ topology::Field<topology::SubMesh>* rateTimeField = (_dbRate) ? &_parameters->get("rate time") : 0;
+ topology::VecVisitorMesh* rateTimeVisitor = (rateTimeField) ? new topology::VecVisitorMesh(*rateTimeField) : 0;
+ PetscScalar* rateTimeArray = (rateTimeVisitor) ? rateTimeVisitor->localArray() : NULL;
- err = PetscSectionGetDof(valuesSection, v, &vdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(valuesSection, v, &voff);CHECK_PETSC_ERROR(err);
- assert(vdof == spaceDim);
- for(PetscInt d = 0; d < vdof; ++d)
- valuesArray[voff+d] = 0.0;
+ topology::Field<topology::SubMesh>* changeField = (_dbChange) ? &_parameters->get("change") : 0;
+ topology::VecVisitorMesh* changeVisitor = (changeField) ? new topology::VecVisitorMesh(*changeField) : 0;
+ PetscScalar* changeArray = (changeVisitor) ? changeVisitor->localArray() : NULL;
+ topology::Field<topology::SubMesh>* changeTimeField = (_dbChange) ? &_parameters->get("change time") : 0;
+ topology::VecVisitorMesh* changeTimeVisitor = (changeTimeField) ? new topology::VecVisitorMesh(*changeTimeField) : 0;
+ PetscScalar* changeTimeArray = (changeTimeVisitor) ? changeTimeVisitor->localArray() : NULL;
+
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ const PetscInt voff = valueVisitor.sectionOffset(v);
+ assert(spaceDim == valueVisitor.sectionDof(v));
+ for (PetscInt d = 0; d < spaceDim; ++d) {
+ valueArray[voff+d] = 0.0;
+ } // for
+
// Contribution from initial value
if (_dbInitial) {
- PetscInt idof, ioff;
-
- err = PetscSectionGetDof(initialSection, v, &idof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(initialSection, v, &ioff);CHECK_PETSC_ERROR(err);
- assert(idof == vdof);
- for(PetscInt d = 0; d < idof; ++d)
- valuesArray[voff+d] += initialArray[ioff+d];
+ assert(initialVisitor);
+ const PetscInt ioff = initialVisitor->sectionOffset(v);
+ assert(spaceDim == initialVisitor->sectionDof(v));
+ for (PetscInt d = 0; d < spaceDim; ++d) {
+ valueArray[voff+d] += initialArray[ioff+d];
+ } // for
} // if
// Contribution from rate of change of value
if (_dbRate) {
- PetscInt rdof, roff, rtdof, rtoff;
+ assert(rateVisitor);
+ const PetscInt roff = rateVisitor->sectionOffset(v);
+ assert(spaceDim == rateVisitor->sectionDof(v));
+ assert(rateTimeVisitor);
+ const PetscInt rtoff = rateTimeVisitor->sectionOffset(v);
+ assert(1 == rateTimeVisitor->sectionDof(v));
- err = PetscSectionGetDof(rateSection, v, &rdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(rateSection, v, &roff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(rateTimeSection, v, &rtdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(rateTimeSection, v, &rtoff);CHECK_PETSC_ERROR(err);
- assert(rdof == vdof);
- assert(rtdof == 1);
-
const PylithScalar tRel = t - rateTimeArray[rtoff];
if (tRel > 0.0) // rate of change integrated over time
- for(PetscInt d = 0; d < spaceDim; ++d)
- valuesArray[voff+d] += rateArray[roff+d] * tRel;
+ for(int iDim = 0; iDim < spaceDim; ++iDim) {
+ valueArray[voff+iDim] += rateArray[roff+iDim] * tRel;
+ } // for
} // if
-
+
// Contribution from change of value
if (_dbChange) {
- PetscInt cdof, coff, ctdof, ctoff;
+ assert(changeVisitor);
+ const PetscInt coff = changeVisitor->sectionOffset(v);
+ assert(spaceDim == changeVisitor->sectionDof(v));
+ const PetscInt ctoff = changeTimeVisitor->sectionOffset(v);
+ assert(1 == changeTimeVisitor->sectionDof(v));
- err = PetscSectionGetDof(changeSection, v, &cdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(changeSection, v, &coff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(changeTimeSection, v, &ctdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(changeTimeSection, v, &ctoff);CHECK_PETSC_ERROR(err);
- assert(cdof == vdof);
- assert(ctdof == 1);
-
const PylithScalar tRel = t - changeTimeArray[ctoff];
if (tRel >= 0) { // change in value over time
- PylithScalar scale = 1.0;
- if (0 != _dbTimeHistory) {
- PylithScalar tDim = tRel * timeScale;
- /* _getNormalizer().dimensionalize(&tDim, 1, timeScale);*/
- const int err = _dbTimeHistory->query(&scale, tDim);
- if (0 != err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << tDim
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- } // if
- for(PetscInt d = 0; d < spaceDim; ++d)
- valuesArray[voff+d] += changeArray[coff+d] * scale;
+ PylithScalar scale = 1.0;
+ if (_dbTimeHistory) {
+ PylithScalar tDim = tRel*_timeScale;
+ const int err = _dbTimeHistory->query(&scale, tDim);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Error querying for time '" << tDim
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+ for (int iDim = 0; iDim < spaceDim; ++iDim) {
+ valueArray[voff+iDim] += changeArray[coff+iDim]*scale;
+ } // for
} // if
} // if
} // for
- err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
- if (_dbInitial) {
- err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
- }
- if (_dbRate) {
- err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
- }
- if (_dbChange) {
- err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
- }
+
+ delete initialVisitor; initialVisitor = 0;
+ delete rateVisitor; rateVisitor = 0;
+ delete rateTimeVisitor; rateTimeVisitor = 0;
+ delete changeVisitor; changeVisitor = 0;
+ delete changeTimeVisitor; changeTimeVisitor = 0;
} // calculate
@@ -454,75 +432,55 @@
assert(_parameters);
// Get vertices.
- DM dmMesh = _parameters->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);
-
- const spatialdata::geocoords::CoordSys* cs = _parameters->mesh().coordsys();
- assert(cs);
+ const spatialdata::geocoords::CoordSys* cs = _parameters->mesh().coordsys();assert(cs);
const int spaceDim = cs->spaceDim();
const PylithScalar lengthScale = normalizer.lengthScale();
- // Containers for database query results and quadrature coordinates in
- // reference geometry.
- scalar_array valuesVertex(querySize);
+ scalar_array coordsVertex(spaceDim);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar *coordArray = coordsVisitor.localArray();
- // Get sections.
- scalar_array coordsVertexGlobal(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::Field<topology::SubMesh>& parametersField = _parameters->get(name);
+ topology::VecVisitorMesh parametersVisitor(parametersField);
+ PetscScalar* parametersArray = parametersVisitor.localArray();
- PetscSection valueSection = _parameters->get(name).petscSection();
- Vec valueVec = _parameters->get(name).localVector();
- PetscScalar *tractionsArray;
- assert(valueSection);assert(valueVec);
-
- // Loop over cells in boundary mesh and perform queries.
- err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(valueVec, &tractionsArray);CHECK_PETSC_ERROR(err);
+ // Containers for database query results and quadrature coordinates in
+ // reference geometry.
+ scalar_array valueVertex(querySize);
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt cdof, coff;
+ // Get dimensionalized coordinates of vertex
+ const int coff = coordsVisitor.sectionOffset(v);
+ assert(spaceDim == coordsVisitor.sectionDof(v));
+ for (PetscInt d = 0; d < spaceDim; ++d) {
+ coordsVertex[d] = coordArray[coff+d];
+ } // for
+ normalizer.dimensionalize(&coordsVertex[0], coordsVertex.size(), lengthScale);
- err = PetscSectionGetDof(coordSection, v, &cdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(coordSection, v, &coff);CHECK_PETSC_ERROR(err);
- assert(spaceDim == cdof);
- for(PetscInt d = 0; d < cdof; ++d) {
- coordsVertexGlobal[d] = coordArray[coff+d];
- }
- normalizer.dimensionalize(&coordsVertexGlobal[0], coordsVertexGlobal.size(), lengthScale);
-
- valuesVertex = 0.0;
- err = db->query(&valuesVertex[0], querySize, &coordsVertexGlobal[0], spaceDim, cs);
+ valueVertex = 0.0;
+ int err = db->query(&valueVertex[0], valueVertex.size(), &coordsVertex[0], coordsVertex.size(), cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find values at (";
+ msg << "Error querying for '" << name << "' at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << coordsVertexGlobal[i];
- msg << ") for traction boundary condition " << _label << "\n"
- << "using spatial database " << db->label() << ".";
+ msg << " " << coordsVertex[i];
+ msg << ") using spatial database " << db->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(&valueVertex[0], valueVertex.size(), scale);
- normalizer.nondimensionalize(&valuesVertex[0], valuesVertex.size(), scale);
-
// Update section
- PetscInt vdof, voff;
-
- err = PetscSectionGetDof(valueSection, v, &vdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(valueSection, v, &voff);CHECK_PETSC_ERROR(err);
- assert(vdof == querySize);
- for(PetscInt d = 0; d < vdof; ++d)
- tractionsArray[voff+d] = valuesVertex[d];
+ const PetscInt off = parametersVisitor.sectionOffset(v);
+ assert(querySize == parametersVisitor.sectionDof(v));
+ for(int i = 0; i < querySize; ++i) {
+ parametersArray[off+i] = valueVertex[i];
+ } // for
} // for
- err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
} // _queryDB
// End of file
More information about the CIG-COMMITS
mailing list