[cig-commits] r20935 - in short/3D/PyLith/trunk: libsrc/pylith/faults unittests/libtests/faults
knepley at geodynamics.org
knepley at geodynamics.org
Thu Oct 25 19:56:16 PDT 2012
Author: knepley
Date: 2012-10-25 19:56:15 -0700 (Thu, 25 Oct 2012)
New Revision: 20935
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/LiuCosSlipFn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/TimeHistorySlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.cc
Log:
More working fault tests
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/LiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/LiuCosSlipFn.cc 2012-10-26 02:13:20 UTC (rev 20934)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/LiuCosSlipFn.cc 2012-10-26 02:56:15 UTC (rev 20935)
@@ -87,27 +87,28 @@
const PylithScalar timeScale = normalizer.timeScale();
// 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();
+ DM dmMesh = faultMesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Fault");
delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
assert(0 != _parameters);
_parameters->add("final slip", "final_slip");
- topology::Field<topology::SubMesh>& finalSlip =
- _parameters->get("final slip");
- finalSlip.newSection(vertices, spaceDim);
+ 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);
- const ALE::Obj<RealSection>& finalSlipSection = finalSlip.section();
- assert(!finalSlipSection.isNull());
+ PetscSection finalSlipSection = finalSlip.petscSection();
+ Vec finalSlipVec = finalSlip.localVector();
+ PetscScalar *finalSlipArray;
+ assert(finalSlipSection);assert(finalSlipVec);
_parameters->add("slip time", "slip_time");
topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
@@ -115,16 +116,21 @@
slipTime.allocate();
slipTime.scale(timeScale);
slipTime.vectorFieldType(topology::FieldBase::SCALAR);
- const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
- assert(!slipTimeSection.isNull());
+ PetscSection slipTimeSection = slipTime.petscSection();
+ Vec slipTimeVec = slipTime.localVector();
+ PetscScalar *slipTimeArray;
+ assert(slipTimeSection);assert(slipTimeVec);
_parameters->add("rise time", "rise_time");
topology::Field<topology::SubMesh>& riseTime = _parameters->get("rise time");
- riseTime.cloneSection(slipTime);
+ riseTime.newSection(finalSlip, 1);
+ riseTime.allocate();
riseTime.scale(timeScale);
riseTime.vectorFieldType(topology::FieldBase::SCALAR);
- const ALE::Obj<RealSection>& riseTimeSection = riseTime.section();
- assert(!riseTimeSection.isNull());
+ PetscSection riseTimeSection = riseTime.petscSection();
+ Vec riseTimeVec = riseTime.localVector();
+ PetscScalar *riseTimeArray;
+ assert(riseTimeSection);assert(riseTimeVec);
logger.stagePop();
@@ -163,40 +169,54 @@
_dbRiseTime->queryVals(riseTimeValues, 1);
// Get coordinates of vertices
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
+ scalar_array vCoordsGlobal(spaceDim);
+ PetscSection coordSection;
+ Vec coordVec;
+ PetscScalar *coordArray;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
_slipVertex.resize(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);
+ 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];
+ }
+ normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(), lengthScale);
- int err = _dbFinalSlip->query(&_slipVertex[0], _slipVertex.size(),
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+ int err = _dbFinalSlip->query(&_slipVertex[0], _slipVertex.size(), &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
msg << "Could not find slip rate at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoordsGlobal[i];
+ msg << " " << vCoordsGlobal[i];
msg << ") using spatial database " << _dbFinalSlip->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(&_slipVertex[0], _slipVertex.size(),
- lengthScale);
+ normalizer.nondimensionalize(&_slipVertex[0], _slipVertex.size(), lengthScale);
- err = _dbSlipTime->query(&_slipTimeVertex, 1,
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+ err = _dbSlipTime->query(&_slipTimeVertex, 1, &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
msg << "Could not find slip initiation time at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoordsGlobal[i];
+ msg << " " << vCoordsGlobal[i];
msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
@@ -204,22 +224,27 @@
// add origin time to rupture time
_slipTimeVertex += originTime;
- err = _dbRiseTime->query(&_riseTimeVertex, 1,
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+ err = _dbRiseTime->query(&_riseTimeVertex, 1, &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
msg << "Could not find rise time at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoordsGlobal[i];
+ msg << " " << vCoordsGlobal[i];
msg << ") using spatial database " << _dbRiseTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
normalizer.nondimensionalize(&_riseTimeVertex, 1, timeScale);
- finalSlipSection->updatePoint(*v_iter, &_slipVertex[0]);
- slipTimeSection->updatePoint(*v_iter, &_slipTimeVertex);
- riseTimeSection->updatePoint(*v_iter, &_riseTimeVertex);
+ for(PetscInt d = 0; d < fsdof; ++d) {
+ finalSlipArray[fsoff+d] = _slipVertex[d];
+ }
+ 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();
@@ -237,53 +262,69 @@
assert(0 != _parameters);
// Get vertices in fault mesh
- const ALE::Obj<SieveMesh>& sieveMesh = slip->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();
+ DM dmMesh = slip->mesh().dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Get sections
- const topology::Field<topology::SubMesh>& finalSlip =
- _parameters->get("final slip");
- const ALE::Obj<RealSection>& finalSlipSection = finalSlip.section();
- assert(!finalSlipSection.isNull());
- const topology::Field<topology::SubMesh>& slipTime =
- _parameters->get("slip time");
- const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
- assert(!slipTimeSection.isNull());
- const topology::Field<topology::SubMesh>& riseTime =
- _parameters->get("rise time");
- const ALE::Obj<RealSection>& riseTimeSection = riseTime.section();
- assert(!riseTimeSection.isNull());
- const ALE::Obj<RealSection>& slipSection = slip->section();
- assert(!slipSection.isNull());
+ const topology::Field<topology::SubMesh>& finalSlip = _parameters->get("final slip");
+ PetscSection finalSlipSection = finalSlip.petscSection();
+ Vec finalSlipVec = finalSlip.localVector();
+ PetscScalar *finalSlipArray;
+ assert(finalSlipSection);assert(finalSlipVec);
+ const topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
+ PetscSection slipTimeSection = slipTime.petscSection();
+ Vec slipTimeVec = slipTime.localVector();
+ PetscScalar *slipTimeArray;
+ assert(slipTimeSection);assert(slipTimeVec);
+ 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);
const int spaceDim = _slipVertex.size();
- for (label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- finalSlipSection->restrictPoint(*v_iter, &_slipVertex[0],
- _slipVertex.size());
- slipTimeSection->restrictPoint(*v_iter, &_slipTimeVertex, 1);
- riseTimeSection->restrictPoint(*v_iter, &_riseTimeVertex, 1);
+ 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;
+ 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);
+
PylithScalar finalSlipMag = 0.0;
for (int i=0; i < spaceDim; ++i)
- finalSlipMag += _slipVertex[i]*_slipVertex[i];
+ finalSlipMag += finalSlipArray[fsoff+i]*finalSlipArray[fsoff+i];
finalSlipMag = sqrt(finalSlipMag);
- const PylithScalar slip = _slipFn(t-_slipTimeVertex, finalSlipMag,
- _riseTimeVertex);
+ const PylithScalar slip = _slipFn(t-slipTimeArray[stoff], finalSlipMag, riseTimeArray[rtoff]);
const PylithScalar scale = finalSlipMag > 0.0 ? slip / finalSlipMag : 0.0;
- _slipVertex *= scale;
// Update field
- slipSection->updateAddPoint(*v_iter, &_slipVertex[0]);
+ for(PetscInt d = 0; d < sdof; ++d) {
+ slipArray[soff+d] += finalSlipArray[fsoff+d] * scale;
+ }
} // for
- PetscLogFlops(vertices->size() * (2+28 + 3*_slipVertex.size()));
+ PetscLogFlops((vEnd-vStart) * (2+28 + 3*_slipVertex.size()));
} // slip
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/TimeHistorySlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/TimeHistorySlipFn.cc 2012-10-26 02:13:20 UTC (rev 20934)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/TimeHistorySlipFn.cc 2012-10-26 02:56:15 UTC (rev 20935)
@@ -93,25 +93,26 @@
logger.stagePush("Fault");
// 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();
+ DM dmMesh = faultMesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
assert(0 != _parameters);
_parameters->add("slip amplitude", "slip_amplitude");
- topology::Field<topology::SubMesh>& slipAmplitude =
- _parameters->get("slip amplitude");
- slipAmplitude.newSection(vertices, spaceDim);
+ topology::Field<topology::SubMesh>& slipAmplitude = _parameters->get("slip amplitude");
+ slipAmplitude.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
slipAmplitude.allocate();
slipAmplitude.scale(lengthScale);
slipAmplitude.vectorFieldType(topology::FieldBase::VECTOR);
- const ALE::Obj<RealSection>& slipAmplitudeSection = slipAmplitude.section();
- assert(!slipAmplitudeSection.isNull());
+ PetscSection finalSlipSection = slipAmplitude.petscSection();
+ Vec finalSlipVec = slipAmplitude.localVector();
+ PetscScalar *finalSlipArray;
+ assert(finalSlipSection);assert(finalSlipVec);
_parameters->add("slip time", "slip_time");
topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
@@ -119,8 +120,10 @@
slipTime.allocate();
slipTime.scale(timeScale);
slipTime.vectorFieldType(topology::FieldBase::SCALAR);
- const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
- assert(!slipTimeSection.isNull());
+ PetscSection slipTimeSection = slipTime.petscSection();
+ Vec slipTimeVec = slipTime.localVector();
+ PetscScalar *slipTimeArray;
+ assert(slipTimeSection);assert(slipTimeVec);
logger.stagePop();
@@ -155,32 +158,44 @@
_dbSlipTime->queryVals(slipTimeValues, 1);
// Get coordinates of vertices
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
+ scalar_array vCoordsGlobal(spaceDim);
+ PetscSection coordSection;
+ Vec coordVec;
+ PetscScalar *coordArray;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
_slipVertex.resize(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);
+ 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];
+ }
+ normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(), lengthScale);
- int err = _dbAmplitude->query(&_slipVertex[0], _slipVertex.size(),
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+ int err = _dbAmplitude->query(&_slipVertex[0], _slipVertex.size(), &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
msg << "Could not find slip amplitude at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoordsGlobal[i];
+ msg << " " << vCoordsGlobal[i];
msg << ") using spatial database " << _dbAmplitude->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(&_slipVertex[0], _slipVertex.size(),
- lengthScale);
+ normalizer.nondimensionalize(&_slipVertex[0], _slipVertex.size(), lengthScale);
err = _dbSlipTime->query(&_slipTimeVertex, 1,
&vCoordsGlobal[0], vCoordsGlobal.size(), cs);
@@ -188,7 +203,7 @@
std::ostringstream msg;
msg << "Could not find slip initiation time at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoordsGlobal[i];
+ msg << " " << vCoordsGlobal[i];
msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
@@ -196,8 +211,10 @@
// add origin time to rupture time
_slipTimeVertex += originTime;
- slipAmplitudeSection->updatePoint(*v_iter, &_slipVertex[0]);
- slipTimeSection->updatePoint(*v_iter, &_slipTimeVertex);
+ for(PetscInt d = 0; d < fsdof; ++d) {
+ finalSlipArray[fsoff+d] = _slipVertex[d];
+ }
+ slipTimeArray[stoff] = _slipTimeVertex;
} // for
// Close databases.
@@ -220,54 +237,64 @@
assert(0 != _dbTimeHistory);
// Get vertices in fault mesh
- const ALE::Obj<SieveMesh>& sieveMesh = slip->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();
+ DM dmMesh = slip->mesh().dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Get sections
- const topology::Field<topology::SubMesh>& slipAmplitude =
- _parameters->get("slip amplitude");
- const ALE::Obj<RealSection>& slipAmplitudeSection = slipAmplitude.section();
- assert(!slipAmplitudeSection.isNull());
- const topology::Field<topology::SubMesh>& slipTime =
- _parameters->get("slip time");
- const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
- assert(!slipTimeSection.isNull());
- const ALE::Obj<RealSection>& slipSection = slip->section();
- assert(!slipSection.isNull());
+ 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>& 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);
PylithScalar amplitude = 0.0;
- for (label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- slipAmplitudeSection->restrictPoint(*v_iter,
- &_slipVertex[0], _slipVertex.size());
- slipTimeSection->restrictPoint(*v_iter, &_slipTimeVertex, 1);
+ 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;
- PylithScalar relTime = t - _slipTimeVertex;
- if (relTime < 0.0)
- _slipVertex = 0.0;
- else {
+ 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);
+
+ PylithScalar relTime = t - slipTimeArray[stoff];
+ if (relTime >= 0.0) {
relTime *= _timeScale;
const int err = _dbTimeHistory->query(&litude, relTime);
if (0 != err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << relTime
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
+ std::ostringstream msg;
+ msg << "Error querying for time '" << relTime
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
} // if
- _slipVertex *= amplitude;
+
+ for(PetscInt d = 0; d < sdof; ++d) {
+ slipArray[soff+d] += finalSlipArray[fsoff+d] * amplitude;
+ }
} // else
-
- // Update field
- slipSection->updateAddPoint(*v_iter, &_slipVertex[0]);
} // for
- PetscLogFlops(vertices->size() * 3);
+ PetscLogFlops((vEnd-vStart) * 3);
} // slip
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc 2012-10-26 02:13:20 UTC (rev 20934)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc 2012-10-26 02:56:15 UTC (rev 20935)
@@ -102,13 +102,14 @@
CPPUNIT_ASSERT(0 != cs);
const int spaceDim = cs->spaceDim();
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = faultMesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
topology::Field<topology::SubMesh> slip(faultMesh);
- slip.newSection(vertices, spaceDim);
+ slip.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
slip.allocate();
const PylithScalar t = 2.134;
@@ -116,12 +117,12 @@
const PylithScalar tolerance = 1.0e-06;
int iPoint = 0;
-
- const ALE::Obj<RealSection>& slipSection = slip.section();
- CPPUNIT_ASSERT(!slipSection.isNull());
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != verticesEnd;
- ++v_iter, ++iPoint) {
+ PetscSection slipSection = slip.petscSection();
+ Vec slipVec = slip.localVector();
+ PetscScalar *slipArray;
+ CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
+ err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
PylithScalar slipMag = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim)
slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
@@ -130,15 +131,15 @@
const PylithScalar tau = slipMag / (exp(1.0) * peakRate);
const PylithScalar t0 = slipTimeE[iPoint];
const PylithScalar slipNorm = 1.0 - exp(-(t-t0)/tau) * (1.0 + (t-t0)/tau);
+ PetscInt dof, off;
- const int fiberDim = slipSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const PylithScalar* vals = slipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
-
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const PylithScalar slipE = finalSlipE[iPoint*spaceDim+iDim] * slipNorm;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, vals[iDim], tolerance);
+ err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
+
+ for(PetscInt d = 0; d < dof; ++d) {
+ const PylithScalar slipE = finalSlipE[iPoint*spaceDim+d] * slipNorm;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, slipArray[off+d], tolerance);
} // for
} // for
} // testSlip
@@ -172,7 +173,8 @@
// Set up coordinates
spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(mesh->dimension());
+ const int spaceDim = mesh->dimension();
+ cs.setSpaceDim(spaceDim);
cs.initialize();
mesh->coordsys(&cs);
@@ -198,9 +200,47 @@
// using create() instead of createParallel().
const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- faultSieveMesh->setRealSection("coordinates",
- sieveMesh->getRealSection("coordinates"));
+ const ALE::Obj<RealSection>& oldCoordSection = sieveMesh->getRealSection("coordinates");
+ faultSieveMesh->setRealSection("coordinates", oldCoordSection);
+ DM dmMesh = faultMesh->dmMesh();
+ IS subpointMap;
+ const PetscInt *points;
+ PetscSection coordSection;
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(dmMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(coordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+ }
+ err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt coordSize;
+
+ err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(mesh->comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(coordSection, v, &off);CHECK_PETSC_ERROR(err);
+ const PetscScalar *oldCoords = oldCoordSection->restrictPoint(points[v]);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coords[off+d] = oldCoords[d];
+ }
+ }
+ err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(dmMesh, coordVec);CHECK_PETSC_ERROR(err);
+
// Setup databases
spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc 2012-10-26 02:13:20 UTC (rev 20934)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc 2012-10-26 02:56:15 UTC (rev 20935)
@@ -241,13 +241,14 @@
CPPUNIT_ASSERT(0 != cs);
const int spaceDim = cs->spaceDim();
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = faultMesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
topology::Field<topology::SubMesh> slip(faultMesh);
- slip.newSection(vertices, spaceDim);
+ slip.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
slip.allocate();
const PylithScalar t = 2.134;
@@ -255,11 +256,12 @@
const PylithScalar tolerance = 1.0e-06;
int iPoint = 0;
- const ALE::Obj<RealSection>& slipSection = slip.section();
- CPPUNIT_ASSERT(!slipSection.isNull());
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != verticesEnd;
- ++v_iter, ++iPoint) {
+ PetscSection slipSection = slip.petscSection();
+ Vec slipVec = slip.localVector();
+ PetscScalar *slipArray;
+ CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
+ err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
PylithScalar slipMag = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim)
slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
@@ -268,15 +270,15 @@
const PylithScalar slipNorm = (slipMag > 0.0) ?
_slipFn(t - slipTimeE[iPoint], slipMag, riseTimeE[iPoint]) / slipMag :
0.0;
+ PetscInt dof, off;
- const int fiberDim = slipSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const PylithScalar* vals = slipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
-
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const PylithScalar slipE = finalSlipE[iPoint*spaceDim+iDim] * slipNorm;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, vals[iDim], tolerance);
+ err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
+
+ for(PetscInt d = 0; d < dof; ++d) {
+ const PylithScalar slipE = finalSlipE[iPoint*spaceDim+d] * slipNorm;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, slipArray[off+d], tolerance);
} // for
} // for
} // testSlip
@@ -329,7 +331,8 @@
// Set up coordinates
spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(mesh->dimension());
+ const int spaceDim = mesh->dimension();
+ cs.setSpaceDim(spaceDim);
cs.initialize();
mesh->coordsys(&cs);
@@ -354,9 +357,47 @@
// using create() instead of createParallel().
const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- faultSieveMesh->setRealSection("coordinates",
- sieveMesh->getRealSection("coordinates"));
+ const ALE::Obj<RealSection>& oldCoordSection = sieveMesh->getRealSection("coordinates");
+ faultSieveMesh->setRealSection("coordinates", oldCoordSection);
+ DM dmMesh = faultMesh->dmMesh();
+ IS subpointMap;
+ const PetscInt *points;
+ PetscSection coordSection;
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(dmMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(coordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+ }
+ err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt coordSize;
+
+ err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(mesh->comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(coordSection, v, &off);CHECK_PETSC_ERROR(err);
+ const PetscScalar *oldCoords = oldCoordSection->restrictPoint(points[v]);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coords[off+d] = oldCoords[d];
+ }
+ }
+ err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(dmMesh, coordVec);CHECK_PETSC_ERROR(err);
+
// Setup databases
spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -425,9 +466,47 @@
// using create() instead of createParallel().
const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- faultSieveMesh->setRealSection("coordinates",
- sieveMesh->getRealSection("coordinates"));
+ const ALE::Obj<RealSection>& oldCoordSection = sieveMesh->getRealSection("coordinates");
+ faultSieveMesh->setRealSection("coordinates", oldCoordSection);
+ DM dmMesh = faultMesh.dmMesh();
+ IS subpointMap;
+ const PetscInt *points;
+ PetscSection coordSection;
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(dmMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(coordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+ }
+ err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt coordSize;
+
+ err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(mesh.comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(coordSection, v, &off);CHECK_PETSC_ERROR(err);
+ const PetscScalar *oldCoords = oldCoordSection->restrictPoint(points[v]);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coords[off+d] = oldCoords[d];
+ }
+ }
+ err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(dmMesh, coordVec);CHECK_PETSC_ERROR(err);
+
// Setup databases
spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -455,46 +534,47 @@
slipfn.initialize(faultMesh, normalizer, originTime);
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
-
CPPUNIT_ASSERT(0 != slipfn._parameters);
- const ALE::Obj<RealSection>& finalSlipSection =
- slipfn._parameters->get("final slip").section();
- CPPUNIT_ASSERT(!finalSlipSection.isNull());
- const ALE::Obj<RealSection>& slipTimeSection =
- slipfn._parameters->get("slip time").section();
- CPPUNIT_ASSERT(!slipTimeSection.isNull());
- const ALE::Obj<RealSection>& riseTimeSection =
- slipfn._parameters->get("rise time").section();
- CPPUNIT_ASSERT(!riseTimeSection.isNull());
+ PetscSection finalSlipSection = slipfn._parameters->get("final slip").petscSection();
+ Vec finalSlipVec = slipfn._parameters->get("final slip").localVector();
+ PetscScalar *finalSlipArray;
+ CPPUNIT_ASSERT(finalSlipSection);CPPUNIT_ASSERT(finalSlipVec);
+ PetscSection slipTimeSection = slipfn._parameters->get("slip time").petscSection();
+ Vec slipTimeVec = slipfn._parameters->get("slip time").localVector();
+ PetscScalar *slipTimeArray;
+ CPPUNIT_ASSERT(slipTimeSection);CPPUNIT_ASSERT(slipTimeVec);
+ PetscSection riseTimeSection = slipfn._parameters->get("rise time").petscSection();
+ Vec riseTimeVec = slipfn._parameters->get("rise time").localVector();
+ PetscScalar *riseTimeArray;
+ CPPUNIT_ASSERT(riseTimeSection);CPPUNIT_ASSERT(riseTimeVec);
const PylithScalar tolerance = 1.0e-06;
int iPoint = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != verticesEnd;
- ++v_iter, ++iPoint) {
- CPPUNIT_ASSERT_EQUAL(spaceDim, finalSlipSection->getFiberDimension(*v_iter));
- const PylithScalar* finalSlipVertex = finalSlipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != finalSlipVertex);
- for (int iDim=0; iDim < spaceDim; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+iDim],
- finalSlipVertex[iDim],
- tolerance);
+ 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, ++iPoint) {
+ PetscInt fsdof, fsoff, stdof, stoff, rtdof, rtoff;
- CPPUNIT_ASSERT_EQUAL(1, slipTimeSection->getFiberDimension(*v_iter));
- const PylithScalar* slipTimeVertex = slipTimeSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != slipTimeVertex);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime,
- slipTimeVertex[0], tolerance);
+ 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);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, fsdof);
+ for(PetscInt d = 0; d < fsdof; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+d], finalSlipArray[fsoff+d], tolerance);
- CPPUNIT_ASSERT_EQUAL(1, riseTimeSection->getFiberDimension(*v_iter));
- const PylithScalar* riseTimeVertex = riseTimeSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != riseTimeVertex);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.riseTimeE[iPoint],
- riseTimeVertex[0], tolerance);
+ CPPUNIT_ASSERT_EQUAL(1, stdof);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
+
+ CPPUNIT_ASSERT_EQUAL(1, rtdof);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.riseTimeE[iPoint], riseTimeArray[rtoff], tolerance);
} // 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);
} // _testInitialize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.cc 2012-10-26 02:13:20 UTC (rev 20934)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestTimeHistorySlipFn.cc 2012-10-26 02:56:15 UTC (rev 20935)
@@ -235,13 +235,14 @@
CPPUNIT_ASSERT(0 != cs);
const int spaceDim = cs->spaceDim();
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = faultMesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
topology::Field<topology::SubMesh> slip(faultMesh);
- slip.newSection(vertices, spaceDim);
+ slip.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
slip.allocate();
const PylithScalar t = 2.0;
@@ -249,19 +250,20 @@
const PylithScalar tolerance = 1.0e-06;
int iPoint = 0;
- const ALE::Obj<RealSection>& slipSection = slip.section();
- CPPUNIT_ASSERT(!slipSection.isNull());
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != verticesEnd;
- ++v_iter, ++iPoint) {
- const int fiberDim = slipSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const PylithScalar* vals = slipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ PetscSection slipSection = slip.petscSection();
+ Vec slipVec = slip.localVector();
+ PetscScalar *slipArray;
+ CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
+ err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
- for (int iDim=0; iDim < fiberDim; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE[iPoint*spaceDim+iDim],
- vals[iDim], tolerance);
+ for(PetscInt d = 0; d < dof; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE[iPoint*spaceDim+d], slipArray[off+d], tolerance);
} // for
} // testSlip
@@ -291,7 +293,8 @@
// Set up coordinates
spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(mesh->dimension());
+ const int spaceDim = mesh->dimension();
+ cs.setSpaceDim(spaceDim);
cs.initialize();
mesh->coordsys(&cs);
@@ -316,9 +319,47 @@
// using create() instead of createParallel().
const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- faultSieveMesh->setRealSection("coordinates",
- sieveMesh->getRealSection("coordinates"));
+ const ALE::Obj<RealSection>& oldCoordSection = sieveMesh->getRealSection("coordinates");
+ faultSieveMesh->setRealSection("coordinates", oldCoordSection);
+ DM dmMesh = faultMesh->dmMesh();
+ IS subpointMap;
+ const PetscInt *points;
+ PetscSection coordSection;
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(dmMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(coordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+ }
+ err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt coordSize;
+
+ err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(mesh->comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(coordSection, v, &off);CHECK_PETSC_ERROR(err);
+ const PetscScalar *oldCoords = oldCoordSection->restrictPoint(points[v]);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coords[off+d] = oldCoords[d];
+ }
+ }
+ err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(dmMesh, coordVec);CHECK_PETSC_ERROR(err);
+
// Setup databases
spatialdata::spatialdb::SimpleDB dbAmplitude("slip amplitude");
spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -385,9 +426,47 @@
// using create() instead of createParallel().
const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- faultSieveMesh->setRealSection("coordinates",
- sieveMesh->getRealSection("coordinates"));
+ const ALE::Obj<RealSection>& oldCoordSection = sieveMesh->getRealSection("coordinates");
+ faultSieveMesh->setRealSection("coordinates", oldCoordSection);
+ DM dmMesh = faultMesh.dmMesh();
+ IS subpointMap;
+ const PetscInt *points;
+ PetscSection coordSection;
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(dmMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(coordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+ }
+ err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt coordSize;
+
+ err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(mesh.comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(coordSection, v, &off);CHECK_PETSC_ERROR(err);
+ const PetscScalar *oldCoords = oldCoordSection->restrictPoint(points[v]);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coords[off+d] = oldCoords[d];
+ }
+ }
+ err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(dmMesh, coordVec);CHECK_PETSC_ERROR(err);
+
// Setup databases
spatialdata::spatialdb::SimpleDB dbAmplitude("slip amplitude");
spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -413,36 +492,33 @@
slipfn.initialize(faultMesh, normalizer, originTime);
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
-
CPPUNIT_ASSERT(0 != slipfn._parameters);
- const ALE::Obj<RealSection>& finalSlipSection =
- slipfn._parameters->get("slip amplitude").section();
- CPPUNIT_ASSERT(!finalSlipSection.isNull());
- const ALE::Obj<RealSection>& slipTimeSection =
- slipfn._parameters->get("slip time").section();
- CPPUNIT_ASSERT(!slipTimeSection.isNull());
+ PetscSection finalSlipSection = slipfn._parameters->get("slip amplitude").petscSection();
+ Vec finalSlipVec = slipfn._parameters->get("slip amplitude").localVector();
+ PetscScalar *finalSlipArray;
+ CPPUNIT_ASSERT(finalSlipSection);CPPUNIT_ASSERT(finalSlipVec);
+ PetscSection slipTimeSection = slipfn._parameters->get("slip time").petscSection();
+ Vec slipTimeVec = slipfn._parameters->get("slip time").localVector();
+ PetscScalar *slipTimeArray;
+ CPPUNIT_ASSERT(slipTimeSection);CPPUNIT_ASSERT(slipTimeVec);
const PylithScalar tolerance = 1.0e-06;
int iPoint = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != verticesEnd;
- ++v_iter, ++iPoint) {
- CPPUNIT_ASSERT_EQUAL(spaceDim, finalSlipSection->getFiberDimension(*v_iter));
- const PylithScalar* amplitudeVertex = finalSlipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != amplitudeVertex);
- for (int iDim=0; iDim < spaceDim; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.amplitudeE[iPoint*spaceDim+iDim],
- amplitudeVertex[iDim],
- tolerance);
+ err = VecGetArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
+ PetscInt fsdof, fsoff, stdof, stoff;
- CPPUNIT_ASSERT_EQUAL(1, slipTimeSection->getFiberDimension(*v_iter));
- const PylithScalar* slipTimeVertex = slipTimeSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != slipTimeVertex);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime,
- slipTimeVertex[0], tolerance);
+ 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);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, fsdof);
+ for(PetscInt d = 0; d < fsdof; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.amplitudeE[iPoint*spaceDim+d], finalSlipArray[fsoff+d], tolerance);
+
+ CPPUNIT_ASSERT_EQUAL(1, stdof);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
} // for
} // _testInitialize
More information about the CIG-COMMITS
mailing list