[cig-commits] r20930 - in short/3D/PyLith/trunk: libsrc/pylith/faults unittests/libtests/faults
knepley at geodynamics.org
knepley at geodynamics.org
Thu Oct 25 12:35:03 PDT 2012
Author: knepley
Date: 2012-10-25 12:35:03 -0700 (Thu, 25 Oct 2012)
New Revision: 20930
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/ConstRateSlipFn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/StepSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
Log:
Now a fault DM is created (not in parallel), Fixed some initial slip models
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc 2012-10-25 18:26:37 UTC (rev 20929)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/BruneSlipFn.cc 2012-10-25 19:35:03 UTC (rev 20930)
@@ -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,73 @@
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
+ 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(vertices->size() * (2+8 + 3*_slipVertex.size()));
+ PetscLogFlops((vEnd-vStart) * (2+8 + 3*spaceDim));
} // slip
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2012-10-25 18:26:37 UTC (rev 20929)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2012-10-25 19:35:03 UTC (rev 20930)
@@ -48,7 +48,7 @@
ALE::Obj<SieveFlexMesh>& faultBoundary,
const topology::Mesh& mesh,
const ALE::Obj<topology::Mesh::IntSection>& groupField,
- const bool flipFault)
+ const bool flipFault)
{ // createFault
assert(0 != faultMesh);
assert(!groupField.isNull());
@@ -59,23 +59,18 @@
faultMesh->coordsys(mesh);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
assert(!sieveMesh.isNull());
- ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
- faultSieveMesh =
- new SieveSubMesh(mesh.comm(), mesh.dimension()-1, mesh.debug());
-
- const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
assert(!sieve.isNull());
- const ALE::Obj<SieveSubMesh::sieve_type> ifaultSieve =
- new SieveMesh::sieve_type(sieve->comm(), sieve->debug());
+ ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
+ faultSieveMesh = new SieveSubMesh(mesh.comm(), mesh.dimension()-1, mesh.debug());
+ const ALE::Obj<SieveSubMesh::sieve_type> ifaultSieve = new SieveMesh::sieve_type(sieve->comm(), sieve->debug());
assert(!ifaultSieve.isNull());
- ALE::Obj<SieveFlexMesh> fault =
- new SieveFlexMesh(mesh.comm(), mesh.dimension()-1, mesh.debug());
- assert(!fault.isNull());
- ALE::Obj<SieveFlexMesh::sieve_type> faultSieve =
- new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
- assert(!faultSieve.isNull());
+
+ ALE::Obj<SieveFlexMesh> fault = new SieveFlexMesh(mesh.comm(), mesh.dimension()-1, mesh.debug());
+ ALE::Obj<SieveFlexMesh::sieve_type> faultSieve = new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
+ assert(!fault.isNull());assert(!faultSieve.isNull());
const int debug = mesh.debug();
// Create set with vertices on fault
@@ -123,6 +118,35 @@
ALE::ISieveConverter::convertMesh(*fault, *faultSieveMesh, renumbering, false);
renumbering.clear();
+ // Convert fault to a DM
+ {
+ DM dm;
+ IS subpointMap;
+ PetscInt *renum;
+ PetscInt pStart, pEnd;
+ PetscErrorCode err;
+ SieveSubMesh::renumbering_type renumbering;
+
+ ALE::ISieveConverter::convertMesh(*fault, &dm, renumbering, true);
+ // Have to make subpointMap here: renumbering[original] = fault
+ err = DMComplexGetChart(dm, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ assert(renumbering.size() == pEnd-pStart);
+ err = PetscMalloc((pEnd-pStart) * sizeof(PetscInt), &renum);CHECK_PETSC_ERROR(err);
+ for(SieveSubMesh::renumbering_type::const_iterator p_iter = renumbering.begin(); p_iter != renumbering.end(); ++p_iter) {
+ renum[p_iter->second] = p_iter->first;
+#if 0
+ std::cout << "renum["<<p_iter->second<<"]: "<<p_iter->first<<std::endl;
+#endif
+ }
+ for(PetscInt p = 1; p < pEnd-pStart; ++p) {
+ assert(renum[p] > renum[p-1]);
+ }
+ err = ISCreateGeneral(fault->comm(), pEnd-pStart, renum, PETSC_OWN_POINTER, &subpointMap);CHECK_PETSC_ERROR(err);
+ err = DMComplexSetSubpointMap(dm, subpointMap);CHECK_PETSC_ERROR(err);
+ renumbering.clear();
+ faultMesh->setDMMesh(dm);
+ }
+
logger.stagePop();
} // createFault
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/ConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/ConstRateSlipFn.cc 2012-10-25 18:26:37 UTC (rev 20929)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/ConstRateSlipFn.cc 2012-10-25 19:35:03 UTC (rev 20930)
@@ -85,13 +85,13 @@
normalizer.lengthScale() / 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");
@@ -99,12 +99,14 @@
assert(0 != _parameters);
_parameters->add("slip rate", "slip_rate");
topology::Field<topology::SubMesh>& slipRate = _parameters->get("slip rate");
- slipRate.newSection(vertices, spaceDim);
+ slipRate.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
slipRate.allocate();
slipRate.scale(velocityScale);
slipRate.vectorFieldType(topology::FieldBase::VECTOR);
- const ALE::Obj<RealSection>& slipRateSection = slipRate.section();
- assert(!slipRateSection.isNull());
+ PetscSection slipRateSection = slipRate.petscSection();
+ Vec slipRateVec = slipRate.localVector();
+ PetscScalar *slipRateArray;
+ assert(slipRateSection);assert(slipRateVec);
_parameters->add("slip time", "slip_time");
topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
@@ -112,8 +114,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();
@@ -148,40 +152,51 @@
_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);
_slipRateVertex.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(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];
+ }
+ normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(), lengthScale);
- int err = _dbSlipRate->query(&_slipRateVertex[0], _slipRateVertex.size(),
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+ int err = _dbSlipRate->query(&_slipRateVertex[0], _slipRateVertex.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 " << _dbSlipRate->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(&_slipRateVertex[0], _slipRateVertex.size(),
- velocityScale);
+ normalizer.nondimensionalize(&_slipRateVertex[0], _slipRateVertex.size(), velocityScale);
- 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
@@ -189,9 +204,14 @@
// add origin time to rupture time
_slipTimeVertex += originTime;
- slipRateSection->updatePoint(*v_iter, &_slipRateVertex[0]);
- slipTimeSection->updatePoint(*v_iter, &_slipTimeVertex);
+ for(PetscInt d = 0; d < srdof; ++d) {
+ slipRateArray[sroff+d] = _slipRateVertex[d];
+ }
+ 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();
@@ -208,41 +228,56 @@
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>& slipRate =
- _parameters->get("slip rate");
- const ALE::Obj<RealSection>& slipRateSection = slipRate.section();
- assert(!slipRateSection.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>& slipRate = _parameters->get("slip rate");
+ PetscSection slipRateSection = slipRate.petscSection();
+ Vec slipRateVec = slipRate.localVector();
+ PetscScalar *slipRateArray;
+ assert(slipRateSection);assert(slipRateVec);
+ 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);
- for (label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- slipRateSection->restrictPoint(*v_iter, &_slipRateVertex[0],
- _slipRateVertex.size());
- slipTimeSection->restrictPoint(*v_iter, &_slipTimeVertex, 1);
+ err = VecGetArray(slipRateVec, &slipRateArray);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 srdof, sroff, stdof, stoff, sdof, soff;
- const PylithScalar relTime = t - _slipTimeVertex;
- const PylithScalar elapsedTime = (relTime > 0) ? relTime : 0.0;
- _slipRateVertex *= elapsedTime; // Convert slip rate to slip
-
- // Update field
- slipSection->updateAddPoint(*v_iter, &_slipRateVertex[0]);
+ 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 PylithScalar relTime = t - slipTimeArray[stoff];
+ if (relTime > 0.0) {
+ for(PetscInt d = 0; d < sdof; ++d) {
+ slipArray[soff+d] += slipRateArray[sroff+d] * relTime; // Convert slip rate to slip
+ }
+ }
} // 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(vertices->size() * (1 + _slipRateVertex.size()));
+ PetscLogFlops((vEnd-vStart) * (1 + _slipRateVertex.size()));
} // slip
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/StepSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/StepSlipFn.cc 2012-10-25 18:26:37 UTC (rev 20929)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/StepSlipFn.cc 2012-10-25 19:35:03 UTC (rev 20930)
@@ -87,24 +87,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("final slip", "final_slip");
-
topology::Field<topology::SubMesh>& finalSlip = _parameters->get("final slip");
- finalSlip.newSection(vertices, spaceDim);
+ 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");
@@ -112,8 +114,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();
@@ -148,40 +152,51 @@
_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 = _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 final slip 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
@@ -189,9 +204,14 @@
// add origin time to rupture time
_slipTimeVertex += originTime;
- finalSlipSection->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
+ 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();
@@ -208,40 +228,56 @@
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 ALE::Obj<RealSection>& slipSection = slip->section();
- assert(!slipSection.isNull());
+ 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);
- 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);
+ 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 PylithScalar relTime = t - _slipTimeVertex;
- if (relTime < 0.0)
- _slipVertex = 0.0;
-
- // Update field
- slipSection->updateAddPoint(*v_iter, &_slipVertex[0]);
+ 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);
+
+ const PylithScalar relTime = t - slipTimeArray[stoff];
+ if (relTime >= 0.0) {
+ for(PetscInt d = 0; d < sdof; ++d) {
+ slipArray[soff+d] += finalSlipArray[fsoff+d];
+ }
+ }
} // 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(vertices->size() * 1);
+ PetscLogFlops((vEnd-vStart) * 1);
} // slip
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc 2012-10-25 18:26:37 UTC (rev 20929)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc 2012-10-25 19:35:03 UTC (rev 20930)
@@ -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);
@@ -269,14 +271,15 @@
(slipMag > 0.0) ? slipMag / (exp(1.0) * peakRate) : 1.0;
const PylithScalar t0 = slipTimeE[iPoint];
const PylithScalar slipNorm = 1.0 - exp(-(t-t0)/tau) * (1.0 + (t-t0)/tau);
- const int fiberDim = slipSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const PylithScalar* vals = slipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ PetscInt dof, off;
- 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
@@ -333,7 +336,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);
@@ -359,9 +363,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;
@@ -431,9 +473,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;
@@ -461,46 +541,46 @@
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;
+ 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/TestConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc 2012-10-25 18:26:37 UTC (rev 20929)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc 2012-10-25 19:35:03 UTC (rev 20930)
@@ -206,13 +206,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::Field<topology::Mesh>::VERTICES_FIELD, spaceDim);
slip.allocate();
const PylithScalar t = 1.234;
@@ -220,23 +221,24 @@
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) {
const PylithScalar t0 = slipTimeE[iPoint];
- const int fiberDim = slipSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const PylithScalar* vals = slipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ PetscInt dof, off;
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const PylithScalar slipE = (t - slipTimeE[iPoint]) > 0.0 ?
- slipRateE[iPoint*spaceDim+iDim] * (t - slipTimeE[iPoint]) : 0.0;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, vals[iDim], tolerance);
- } // for
+ 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 = (t - t0) > 0.0 ? slipRateE[iPoint*spaceDim+d] * (t - t0) : 0.0;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, slipArray[off+d], tolerance);
+ }
} // for
+ err = VecRestoreArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
} // testSlip
// ----------------------------------------------------------------------
@@ -265,7 +267,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);
@@ -291,9 +294,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 dbSlipRate("slip rate");
spatialdata::spatialdb::SimpleIOAscii ioSlipRate;
@@ -357,9 +398,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 dbSlipRate("slip rate");
spatialdata::spatialdb::SimpleIOAscii ioSlipRate;
@@ -381,37 +460,36 @@
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>& slipRateSection =
- slipfn._parameters->get("slip rate").section();
- CPPUNIT_ASSERT(!slipRateSection.isNull());
- const ALE::Obj<RealSection>& slipTimeSection =
- slipfn._parameters->get("slip time").section();
- CPPUNIT_ASSERT(!slipTimeSection.isNull());
+ PetscSection slipRateSection = slipfn._parameters->get("slip rate").petscSection();
+ Vec slipRateVec = slipfn._parameters->get("slip rate").localVector();
+ PetscScalar *slipRateArray;
+ CPPUNIT_ASSERT(slipRateSection);CPPUNIT_ASSERT(slipRateVec);
+ 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;\
+ 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, slipRateSection->getFiberDimension(*v_iter));
- const PylithScalar* slipRateVertex = slipRateSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != slipRateVertex);
- for (int iDim=0; iDim < spaceDim; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipRateE[iPoint*spaceDim+iDim],
- slipRateVertex[iDim],
- tolerance);
+ err = VecGetArray(slipRateVec, &slipRateArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
+ PetscInt srdof, sroff, 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(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);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, srdof);
+ for(PetscInt d = 0; d < spaceDim; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipRateE[iPoint*spaceDim+d], slipRateArray[sroff+d], tolerance);
+
+ CPPUNIT_ASSERT_EQUAL(1, stdof);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[stoff], tolerance);
} // for
+ err = VecRestoreArray(slipRateVec, &slipRateArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
} // _testInitialize
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc 2012-10-25 18:26:37 UTC (rev 20929)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc 2012-10-25 19:35:03 UTC (rev 20930)
@@ -263,7 +263,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);
@@ -289,9 +290,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;
@@ -351,19 +390,51 @@
data.faultId,
firstFaultVertex, firstLagrangeVertex, firstFaultCell,
useLagrangeConstraints);
- // Need to copy coordinates from mesh to fault mesh since we are not
+ // Need to copy coordinates from mesh to fault mesh since we are
// using create() instead of createParallel().
const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- faultSieveMesh->setRealSection("coordinates",
- sieveMesh->getRealSection("coordinates"));
- DM dmMesh = faultMesh.dmMesh();
- PetscErrorCode err;
+ 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);
- PetscInt vStart, vEnd;
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;
More information about the CIG-COMMITS
mailing list