[cig-commits] r20821 - in short/3D/PyLith/trunk: libsrc/pylith/meshio libsrc/pylith/topology unittests/libtests/bc unittests/libtests/faults
knepley at geodynamics.org
knepley at geodynamics.org
Thu Oct 11 11:29:56 PDT 2012
Author: knepley
Date: 2012-10-11 11:29:55 -0700 (Thu, 11 Oct 2012)
New Revision: 20821
Modified:
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
Log:
Fixes for HDF5Ext
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2012-10-10 17:25:38 UTC (rev 20820)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2012-10-11 18:29:55 UTC (rev 20821)
@@ -248,21 +248,21 @@
PetscInt cStart, cEnd, cMax, vStart, numIndices;
err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- err = DMComplexGetHeightStratum(dmMesh, 0, &vStart, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, PETSC_NULL);CHECK_PETSC_ERROR(err);
err = DMComplexGetVTKBounds(dmMesh, &cMax, PETSC_NULL);CHECK_PETSC_ERROR(err);
if (cMax >= 0) {cEnd = PetscMin(cEnd, cMax);}
err = DMComplexGetCellNumbering(dmMesh, &globalCellNumbers);CHECK_PETSC_ERROR(err);
err = DMComplexGetVertexNumbering(dmMesh, &globalVertexNumbers);CHECK_PETSC_ERROR(err);
- if (cEnd-cStart) {err = DMComplexGetConeSize(dmMesh, cStart, &numCornersLocal);CHECK_PETSC_ERROR(err);}
- err = MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
- err = MPI_Allreduce(&numCellsLocal, &numCells, 1, MPIU_INT, MPI_SUM, ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
err = ISGetSize(globalCellNumbers, &numIndices);CHECK_PETSC_ERROR(err);
err = ISGetIndices(globalCellNumbers, &cells);CHECK_PETSC_ERROR(err);
err = ISGetIndices(globalVertexNumbers, &vertices);CHECK_PETSC_ERROR(err);
for(PetscInt i = 0; i < numIndices; ++i) {
if (cells[i] >= 0) ++numCellsLocal;
}
+ if (cEnd-cStart) {err = DMComplexGetConeSize(dmMesh, cStart, &numCornersLocal);CHECK_PETSC_ERROR(err);}
+ err = MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
+ err = MPI_Allreduce(&numCellsLocal, &numCells, 1, MPIU_INT, MPI_SUM, ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
PetscScalar *tmpVertices = 0;
PetscInt conesSize = numCellsLocal*numCorners;
@@ -369,21 +369,13 @@
try {
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
- const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = field.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- DM dmMesh = mesh.dmMesh();
- assert(dmMesh);
+ DM dmMesh = mesh.dmMesh();
PetscMPIInt commRank;
- PetscErrorCode err = MPI_Comm_rank(((PetscObject) dmMesh)->comm, &commRank);CHECK_PETSC_ERROR(err);
+ PetscErrorCode err;
- const std::string labelName =
- (sieveMesh->hasLabel("censored depth")) ? "censored depth" : "depth";
- ALE::Obj<numbering_type> vNumbering =
- sieveMesh->hasLabel("censored depth") ?
- sieveMesh->getFactory()->getNumbering(sieveMesh, "censored depth", 0) :
- sieveMesh->getFactory()->getNumbering(sieveMesh, 0);
- assert(!vNumbering.isNull());
- field.createScatterWithBC(mesh, vNumbering, context);
+ assert(dmMesh);
+ err = MPI_Comm_rank(((PetscObject) dmMesh)->comm, &commRank);CHECK_PETSC_ERROR(err);
+ field.createScatterWithBC(mesh, PETSC_NULL, context);
field.scatterSectionToVector(context);
PetscVec vector = field.vector(context);
assert(vector);
@@ -399,11 +391,10 @@
binaryViewer = _datasets[field.label()].viewer;
} else {
err = PetscViewerBinaryOpen(((PetscObject) dmMesh)->comm,
- _datasetFilename(field.label()).c_str(),
- FILE_MODE_WRITE, &binaryViewer);
+ _datasetFilename(field.label()).c_str(),
+ FILE_MODE_WRITE, &binaryViewer);
CHECK_PETSC_ERROR(err);
- err = PetscViewerBinarySetSkipHeader(binaryViewer, PETSC_TRUE);
- CHECK_PETSC_ERROR(err);
+ err = PetscViewerBinarySetSkipHeader(binaryViewer, PETSC_TRUE);CHECK_PETSC_ERROR(err);
ExternalDataset dataset;
dataset.numTimeSteps = 0;
dataset.viewer = binaryViewer;
@@ -413,21 +404,30 @@
} // else
assert(binaryViewer);
- err = VecView(vector, binaryViewer);
- CHECK_PETSC_ERROR(err);
+ err = VecView(vector, binaryViewer);CHECK_PETSC_ERROR(err);
++_datasets[field.label()].numTimeSteps;
PetscSection section = field.petscSection();
- PetscInt dof = 0;
+ PetscInt dof = 0, n, numLocalVertices = 0, numVertices;
+ IS globalVertexNumbers;
+
assert(section);
- assert(!sieveMesh->getLabelStratum(labelName, 0).isNull());
- if (sieveMesh->getLabelStratum(labelName, 0)->size() > 0) {
- err = PetscSectionGetDof(section, *sieveMesh->getLabelStratum(labelName, 0)->begin(), &dof);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetVertexNumbering(dmMesh, &globalVertexNumbers);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(globalVertexNumbers, &n);CHECK_PETSC_ERROR(err);
+ if (n > 0) {
+ const PetscInt *indices;
+ err = ISGetIndices(globalVertexNumbers, &indices);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(section, indices[0], &dof);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = 0; v < n; ++v) {
+ if (indices[v] >= 0) ++numLocalVertices;
+ }
+ err = ISRestoreIndices(globalVertexNumbers, &indices);CHECK_PETSC_ERROR(err);
}
int fiberDimLocal = dof;
int fiberDim = 0;
- MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, ((PetscObject) dmMesh)->comm);
- assert(fiberDim > 0);
+ err = MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
+ err = MPI_Allreduce(&numLocalVertices, &numVertices, 1, MPI_INT, MPI_SUM, ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
+ assert(fiberDim > 0);assert(numVertices > 0);
if (!commRank) {
if (createdExternalDataset) {
@@ -437,10 +437,10 @@
hsize_t maxDims[3];
if (3 == ndims) {
maxDims[0] = H5S_UNLIMITED;
- maxDims[1] = vNumbering->getGlobalSize();
+ maxDims[1] = numVertices;
maxDims[2] = fiberDim;
} else {
- maxDims[0] = vNumbering->getGlobalSize();
+ maxDims[0] = numVertices;
maxDims[1] = fiberDim;
} // else
// Create 'vertex_fields' group if necessary.
@@ -463,7 +463,7 @@
const hsize_t ndims = 3;
hsize_t dims[3];
dims[0] = numTimeSteps; // update to current value
- dims[1] = vNumbering->getGlobalSize();
+ dims[1] = numVertices;
dims[2] = fiberDim;
_h5->extendDatasetRawExternal("/vertex_fields", field.label(),
dims, ndims);
@@ -503,26 +503,13 @@
try {
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
- const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = field.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- DM dmMesh = field.mesh().dmMesh();
- assert(dmMesh);
+ DM dmMesh = field.mesh().dmMesh();
PetscMPIInt commRank;
- PetscErrorCode err = MPI_Comm_rank(((PetscObject) dmMesh)->comm, &commRank);CHECK_PETSC_ERROR(err);
+ PetscErrorCode err;
- int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
- int cellDepth = 0;
- err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
- ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
- const int depth = (0 == label) ? cellDepth : labelId;
- const std::string labelName = (0 == label) ?
- ((sieveMesh->hasLabel("censored depth")) ?
- "censored depth" : "depth") : label;
- assert(!sieveMesh->getFactory().isNull());
- const ALE::Obj<typename mesh_type::SieveMesh::numbering_type>& numbering =
- sieveMesh->getFactory()->getNumbering(sieveMesh, labelName, depth);
- assert(!numbering.isNull());
- field.createScatterWithBC(field.mesh(), numbering, context);
+ assert(dmMesh);
+ err = MPI_Comm_rank(((PetscObject) dmMesh)->comm, &commRank);CHECK_PETSC_ERROR(err);
+ field.createScatterWithBC(field.mesh(), PETSC_NULL, context);
field.scatterSectionToVector(context);
PetscVec vector = field.vector(context);
assert(vector);
@@ -538,11 +525,10 @@
binaryViewer = _datasets[field.label()].viewer;
} else {
err = PetscViewerBinaryOpen(((PetscObject) dmMesh)->comm,
- _datasetFilename(field.label()).c_str(),
- FILE_MODE_WRITE, &binaryViewer);
+ _datasetFilename(field.label()).c_str(),
+ FILE_MODE_WRITE, &binaryViewer);
CHECK_PETSC_ERROR(err);
- err = PetscViewerBinarySetSkipHeader(binaryViewer, PETSC_TRUE);
- CHECK_PETSC_ERROR(err);
+ err = PetscViewerBinarySetSkipHeader(binaryViewer, PETSC_TRUE);CHECK_PETSC_ERROR(err);
ExternalDataset dataset;
dataset.numTimeSteps = 0;
dataset.viewer = binaryViewer;
@@ -552,66 +538,75 @@
} // else
assert(binaryViewer);
- err = VecView(vector, binaryViewer);
- CHECK_PETSC_ERROR(err);
+ err = VecView(vector, binaryViewer);CHECK_PETSC_ERROR(err);
++_datasets[field.label()].numTimeSteps;
PetscSection section = field.petscSection();
- PetscInt dof = 0;
+ PetscInt dof = 0, n, numLocalCells = 0, numCells;
+ IS globalCellNumbers;
+
+ /* TODO: Add code for handling a label */
assert(section);
- assert(!sieveMesh->getLabelStratum(labelName, 0).isNull());
- if (sieveMesh->getLabelStratum(labelName, 0)->size() > 0) {
- err = PetscSectionGetDof(section, *sieveMesh->getLabelStratum(labelName, 0)->begin(), &dof);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCellNumbering(dmMesh, &globalCellNumbers);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(globalCellNumbers, &n);CHECK_PETSC_ERROR(err);
+ if (n > 0) {
+ const PetscInt *indices;
+ err = ISGetIndices(globalCellNumbers, &indices);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(section, indices[0], &dof);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = 0; v < n; ++v) {
+ if (indices[v] >= 0) ++numLocalCells;
+ }
+ err = ISRestoreIndices(globalCellNumbers, &indices);CHECK_PETSC_ERROR(err);
}
int fiberDimLocal = dof;
int fiberDim = 0;
MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, ((PetscObject) dmMesh)->comm);
- assert(fiberDim > 0);
+ err = MPI_Allreduce(&numLocalCells, &numCells, 1, MPI_INT, MPI_SUM, ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
+ assert(fiberDim > 0);assert(numCells > 0);
if (!commRank) {
if (createdExternalDataset) {
// Add new external dataset to HDF5 file.
- const int numTimeSteps =
- DataWriter<mesh_type, field_type>::_numTimeSteps;
- const hsize_t ndims = (numTimeSteps > 0) ? 3 : 2;
- hsize_t maxDims[3];
- if (3 == ndims) {
- maxDims[0] = H5S_UNLIMITED;
- maxDims[1] = numbering->getGlobalSize();
- maxDims[2] = fiberDim;
- } else {
- maxDims[0] = numbering->getGlobalSize();
- maxDims[1] = fiberDim;
- } // else
- // Create 'cell_fields' group if necessary.
- if (!_h5->hasGroup("/cell_fields"))
- _h5->createGroup("/cell_fields");
+ const int numTimeSteps =
+ DataWriter<mesh_type, field_type>::_numTimeSteps;
+ const hsize_t ndims = (numTimeSteps > 0) ? 3 : 2;
+ hsize_t maxDims[3];
+ if (3 == ndims) {
+ maxDims[0] = H5S_UNLIMITED;
+ maxDims[1] = numCells;
+ maxDims[2] = fiberDim;
+ } else {
+ maxDims[0] = numCells;
+ maxDims[1] = fiberDim;
+ } // else
+ // Create 'cell_fields' group if necessary.
+ if (!_h5->hasGroup("/cell_fields"))
+ _h5->createGroup("/cell_fields");
- _h5->createDatasetRawExternal("/cell_fields", field.label(),
- _datasetFilename(field.label()).c_str(),
- maxDims, ndims, scalartype);
- std::string fullName = std::string("/cell_fields/") + field.label();
- _h5->writeAttribute(fullName.c_str(), "vector_field_type",
- topology::FieldBase::vectorFieldString(field.vectorFieldType()));
+ _h5->createDatasetRawExternal("/cell_fields", field.label(),
+ _datasetFilename(field.label()).c_str(),
+ maxDims, ndims, scalartype);
+ std::string fullName = std::string("/cell_fields/") + field.label();
+ _h5->writeAttribute(fullName.c_str(), "vector_field_type",
+ topology::FieldBase::vectorFieldString(field.vectorFieldType()));
} else {
- // Update number of time steps in external dataset info in HDF5 file.
- const int totalNumTimeSteps =
- DataWriter<mesh_type, field_type>::_numTimeSteps;
- assert(totalNumTimeSteps > 0);
- const int numTimeSteps = _datasets[field.label()].numTimeSteps;
+ // Update number of time steps in external dataset info in HDF5 file.
+ const int totalNumTimeSteps =
+ DataWriter<mesh_type, field_type>::_numTimeSteps;
+ assert(totalNumTimeSteps > 0);
+ const int numTimeSteps = _datasets[field.label()].numTimeSteps;
- const hsize_t ndims = 3;
- hsize_t dims[3];
- dims[0] = numTimeSteps; // update to current value
- dims[1] = numbering->getGlobalSize();
- dims[2] = fiberDim;
- _h5->extendDatasetRawExternal("/cell_fields", field.label(),
- dims, ndims);
+ const hsize_t ndims = 3;
+ hsize_t dims[3];
+ dims[0] = numTimeSteps; // update to current value
+ dims[1] = numCells;
+ dims[2] = fiberDim;
+ _h5->extendDatasetRawExternal("/cell_fields", field.label(), dims, ndims);
} // if/else
// Update time stamp in "/time, if necessary.
if (_tstampIndex+1 == _datasets[field.label()].numTimeSteps)
- _writeTimeStamp(t);
+ _writeTimeStamp(t);
} // if
} catch (const std::exception& err) {
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-10-10 17:25:38 UTC (rev 20820)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-10-11 18:29:55 UTC (rev 20821)
@@ -1437,7 +1437,6 @@
const typename ALE::Obj<typename SieveMesh::numbering_type> numbering,
const char* context)
{ // createScatterWithBC
- assert(!numbering.isNull());
assert(context);
PetscErrorCode err = 0;
@@ -1454,7 +1453,7 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("GlobalOrder");
- if (!_section.isNull()) {
+ if (!_section.isNull() && !numbering.isNull()) {
// Get global order (create if necessary).
const std::string& orderLabel =
(strlen(context) > 0) ?
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2012-10-10 17:25:38 UTC (rev 20820)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2012-10-11 18:29:55 UTC (rev 20821)
@@ -243,8 +243,8 @@
const PetscInt offset = numCells;
int iConstraint = 0;
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt *cInd, *fcInd;
- PetscInt dof, cdof, fdof, fcdof;
+ const PetscInt *cInd, *fcInd;
+ PetscInt dof, cdof, fdof, fcdof;
err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc 2012-10-10 17:25:38 UTC (rev 20820)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc 2012-10-11 18:29:55 UTC (rev 20821)
@@ -145,8 +145,8 @@
int iConstraint = 0;
for(PetscInt v = vStart; v < vEnd; ++v) {
const int numConstrainedDOF = _data->constraintSizes[v-offset];
- PetscInt *cInd;
- PetscInt dof, cdof;
+ const PetscInt *cInd;
+ PetscInt dof, cdof;
err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc 2012-10-10 17:25:38 UTC (rev 20820)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc 2012-10-11 18:29:55 UTC (rev 20821)
@@ -205,34 +205,36 @@
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;
slipfn.slip(&slip, originTime+t);
const PylithScalar tolerance = 1.0e-06;
+ PetscScalar *slipArray;
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();
+ 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;
- for (int iDim=0; iDim < fiberDim; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE[iPoint*spaceDim+iDim], 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)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE[iPoint*spaceDim+d], slipArray[off+d], tolerance);
} // for
+ err = VecRestoreArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
} // testSlip
// ----------------------------------------------------------------------
@@ -355,7 +357,13 @@
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
faultSieveMesh->setRealSection("coordinates",
sieveMesh->getRealSection("coordinates"));
+ DM dmMesh = faultMesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Setup databases
spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -377,37 +385,35 @@
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());
+ PetscSection finalSlipSection = slipfn._parameters->get("final slip").petscSection();
+ Vec finalSlipVec = slipfn._parameters->get("final slip").localVector();
+ CPPUNIT_ASSERT(finalSlipSection);CPPUNIT_ASSERT(finalSlipVec);
+ PetscSection slipTimeSection = slipfn._parameters->get("slip time").petscSection();
+ Vec slipTimeVec = slipfn._parameters->get("slip time").localVector();
+ CPPUNIT_ASSERT(slipTimeSection);CPPUNIT_ASSERT(slipTimeVec);
const PylithScalar tolerance = 1.0e-06;
+ PetscScalar *finalSlipArray, *slipTimeArray;
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);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
+ PetscInt dof, off;
- 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, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(finalSlipSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
+ for(PetscInt d = 0; d < spaceDim; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+d], finalSlipArray[off+d], tolerance);
+
+ err = PetscSectionGetDof(slipTimeSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(slipTimeSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(1, dof);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime, slipTimeArray[off], tolerance);
} // for
+ err = VecRestoreArray(finalSlipVec, &finalSlipArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(slipTimeVec, &slipTimeArray);CHECK_PETSC_ERROR(err);
} // _testInitialize
More information about the CIG-COMMITS
mailing list