[cig-commits] r22110 - in short/3D/PyLith/trunk/libsrc/pylith: faults feassemble materials meshio
brad at geodynamics.org
brad at geodynamics.org
Mon May 20 10:32:43 PDT 2013
Author: brad
Date: 2013-05-20 10:32:42 -0700 (Mon, 20 May 2013)
New Revision: 22110
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.hh
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/Material.hh
short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilter.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilter.hh
short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.hh
Log:
Added more caching of index sets. Improved efficiency of DataWriterHDF5 write methods (cache number of points and fiber dimension).
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -1685,9 +1685,9 @@
// Get cohesive cells
PetscDM dmMesh = fields.mesh().dmMesh();assert(dmMesh);
- topology::StratumIS cellsStratum(dmMesh, "material-id", id());
- const PetscInt *cellsCohesive = cellsStratum.points();
- const PetscInt numCohesiveCells = cellsStratum.size();
+ assert(_cohesiveIS);
+ const PetscInt *cellsCohesive = _cohesiveIS->points();
+ const PetscInt numCohesiveCells = _cohesiveIS->size();
topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
const PetscInt vStart = verticesStratum.begin();
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -53,7 +53,8 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::faults::FaultCohesiveLagrange::FaultCohesiveLagrange(void)
+pylith::faults::FaultCohesiveLagrange::FaultCohesiveLagrange(void) :
+ _cohesiveIS(0)
{ // constructor
_useLagrangeConstraints = true;
} // constructor
@@ -73,6 +74,7 @@
PYLITH_METHOD_BEGIN;
FaultCohesive::deallocate();
+ delete _cohesiveIS; _cohesiveIS = 0;
PYLITH_METHOD_END;
} // deallocate
@@ -933,9 +935,9 @@
// Check quadrature against mesh
const int numCorners = _quadrature->numBasis();
- topology::StratumIS faultIS(dmMesh, "material-id", id());
- const PetscInt* cells = faultIS.points();
- const PetscInt ncells = faultIS.size();
+ topology::StratumIS cohesiveIS(dmMesh, "material-id", id());
+ const PetscInt* cells = cohesiveIS.points();
+ const PetscInt ncells = cohesiveIS.size();
for(PetscInt i = 0; i < ncells; ++i) {
PetscInt *closure = NULL;
PetscInt cellNumCorners = 0, closureSize;
@@ -1017,9 +1019,9 @@
// Get cohesive cells
PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
- topology::StratumIS faultIS(dmMesh, "material-id", id());
- const PetscInt* cells = faultIS.points();
- const int ncells = faultIS.size();
+ delete _cohesiveIS; _cohesiveIS = new topology::StratumIS(dmMesh, "material-id", id());assert(_cohesiveIS);
+ const PetscInt* cells = _cohesiveIS->points();
+ const int ncells = _cohesiveIS->size();
// Get domain vertices
topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.hh 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.hh 2013-05-20 17:32:42 UTC (rev 22110)
@@ -304,6 +304,8 @@
/// Map label of cohesive cell to label of cells in fault mesh.
std::map<PetscInt, PetscInt> _cohesiveToFault;
+ topology::StratumIS* _cohesiveIS; ///< Index set of cohesive cells.
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -141,8 +141,10 @@
PetscInt coneSize = 0;
for (PetscInt i=0; i < ncells; ++i) {
err = DMPlexGetConeSize(dmMesh, cells[i], &coneSize);PYLITH_CHECK_ERROR(err);
- // TODO: Why isn't this 3? Should be changed to Closure()
+ // TODO: Should be changed to Closure()
if (2*numCorners != coneSize) {
+ // No Lagrange vertices, just negative and positive sides of the
+ // fault, so coneSize is 2*numCorners.
std::ostringstream msg;
msg << "Number of vertices in reference cell (" << numCorners
<< ") is not compatible with number of vertices (" << coneSize
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -105,9 +105,9 @@
// Get cell information
PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
- topology::StratumIS materialIS(dmMesh, "material-id", _material->id());
- const PetscInt* cells = materialIS.points();
- const PetscInt numCells = materialIS.size();
+ assert(_materialIS);
+ const PetscInt* cells = _materialIS->points();
+ const PetscInt numCells = _materialIS->size();
scalar_array dispTCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
@@ -188,9 +188,9 @@
// Get cell information
PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
- topology::StratumIS materialIS(dmMesh, "material-id", _material->id());
- const PetscInt* cells = materialIS.points();
- const PetscInt numCells = materialIS.size();
+ assert(_materialIS);
+ const PetscInt* cells = _materialIS->points();
+ const PetscInt numCells = _materialIS->size();
// Setup field visitors.
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -287,9 +287,9 @@
// Get cells associated with material
PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
- topology::StratumIS materialIS(dmMesh, "material-id", id());
- const PetscInt* cells = materialIS.points();
- const PetscInt numCells = materialIS.size();
+ assert(_materialIS);
+ const PetscInt* cells = _materialIS->points();
+ const PetscInt numCells = _materialIS->size();
// Setup field if necessary.
topology::VecVisitorMesh* fieldVisitor = (field) ? new topology::VecVisitorMesh(*field) : 0;
@@ -376,9 +376,9 @@
// Get cells associated with material
PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
- topology::StratumIS materialIS(dmMesh, "material-id", id());
- const PetscInt* cells = materialIS.points();
- const PetscInt numCells = materialIS.size();
+ assert(_materialIS);
+ const PetscInt* cells = _materialIS->points();
+ const PetscInt numCells = _materialIS->size();
// Setup field if necessary.
topology::VecVisitorMesh* fieldVisitor = (field) ? new topology::VecVisitorMesh(*field) : 0;
@@ -466,27 +466,23 @@
if (field) {
const int numQuadPts = _numQuadPts;
- PetscSection fieldSection = field->petscSection();
- PetscVec fieldVec = field->localVector();
- PetscScalar *fieldArray;
// Get cells associated with material
- PetscDM dmMesh = mesh.dmMesh();
- assert(dmMesh);
- PetscIS cellIS;
- const PetscInt *cells;
- PetscInt numCells;
- PetscErrorCode err = 0;
+ PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
- err = DMPlexGetStratumIS(dmMesh, "material-id", id(), &cellIS);PYLITH_CHECK_ERROR(err);
- err = ISGetLocalSize(cellIS, &numCells);PYLITH_CHECK_ERROR(err);
- err = ISGetIndices(cellIS, &cells);PYLITH_CHECK_ERROR(err);
+ assert(_materialIS);
+ const PetscInt *cells = _materialIS->points();
+ const PetscInt numCells = _materialIS->size();
+ // Setup field if necessary.
+ topology::VecVisitorMesh fieldVisitor(*field);
const int fiberDim = 1*numQuadPts;
bool useCurrentField = false;
- if (fieldSection) {
+ if (fieldVisitor.petscSection()) {
// check fiber dimension
PetscInt fiberDimCurrentLocal = 0;
- if (numCells > 0) {err = PetscSectionGetDof(fieldSection, cells[0], &fiberDimCurrentLocal);PYLITH_CHECK_ERROR(err);}
+ if (numCells > 0) {
+ fiberDimCurrentLocal = fieldVisitor.sectionDof(cells[0]);
+ } // if
PetscInt fiberDimCurrent = 0;
MPI_Allreduce(&fiberDimCurrentLocal, &fiberDimCurrent, 1, MPIU_INT, MPI_MAX, field->mesh().comm());
assert(fiberDimCurrent > 0);
@@ -495,33 +491,24 @@
if (!useCurrentField) {
field->newSection(cells, numCells, fiberDim);
field->allocate();
- fieldSection = field->petscSection();
- fieldVec = field->localVector();
} // if
- assert(fieldSection);
field->label("stable_dt_implicit");
assert(_normalizer);
field->scale(_normalizer->timeScale());
field->vectorFieldType(topology::FieldBase::MULTI_SCALAR);
+ PetscScalar* fieldArray = fieldVisitor.localArray();
scalar_array dtStableCell(numQuadPts);
dtStableCell = PYLITH_MAXSCALAR;
- err = VecGetArray(fieldVec, &fieldArray);PYLITH_CHECK_ERROR(err);
for (PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
- PetscInt dof, off;
- assert(fieldSection);
- err = PetscSectionGetDof(fieldSection, cell, &dof);PYLITH_CHECK_ERROR(err);
- err = PetscSectionGetOffset(fieldSection, cell, &off);PYLITH_CHECK_ERROR(err);
- assert(numQuadPts == dof);
- for (PetscInt d = 0; d < dof; ++d) {
+ const PetscInt off = fieldVisitor.sectionOffset(cell);
+ assert(numQuadPts == fieldVisitor.sectionDof(cell));
+ for (PetscInt d = 0; d < numQuadPts; ++d) {
fieldArray[off+d] = dtStableCell[d];
- }
+ } // for
} // for
- err = VecRestoreArray(fieldVec, &fieldArray);PYLITH_CHECK_ERROR(err);
- err = ISRestoreIndices(cellIS, &cells);PYLITH_CHECK_ERROR(err);
- err = ISDestroy(&cellIS);PYLITH_CHECK_ERROR(err);
} // if
PYLITH_METHOD_RETURN(dtStable);
@@ -578,9 +565,9 @@
// Get cells associated with material
PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
- topology::StratumIS materialIS(dmMesh, "material-id", id());
- const PetscInt* cells = materialIS.points();
- const PetscInt numCells = materialIS.size();
+ assert(_materialIS);
+ const PetscInt* cells = _materialIS->points();
+ const PetscInt numCells = _materialIS->size();
topology::CoordsVisitor coordsVisitor(dmMesh);
@@ -711,9 +698,9 @@
// Get cells associated with material
PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
- topology::StratumIS materialIS(dmMesh, "material-id", id());
- const PetscInt* cells = materialIS.points();
- const PetscInt numCells = materialIS.size();
+ assert(_materialIS);
+ const PetscInt* cells = _materialIS->points();
+ const PetscInt numCells = _materialIS->size();
topology::CoordsVisitor coordsVisitor(dmMesh);
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -45,6 +45,7 @@
_properties(0),
_stateVars(0),
_normalizer(new spatialdata::units::Nondimensional),
+ _materialIS(0),
_numPropsQuadPt(0),
_numVarsQuadPt(0),
_dimension(dimension),
@@ -83,6 +84,7 @@
PYLITH_METHOD_BEGIN;
delete _normalizer; _normalizer = 0;
+ delete _materialIS; _materialIS = 0;
delete _properties; _properties = 0;
delete _stateVars; _stateVars = 0;
@@ -125,9 +127,9 @@
// Get cells associated with material
PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
- topology::StratumIS materialIS(dmMesh, "material-id", _id);
- const PetscInt numCells = materialIS.size();
- const PetscInt* cells = materialIS.points();
+ delete _materialIS; _materialIS = new topology::StratumIS(dmMesh, "material-id", _id);assert(_materialIS);
+ const PetscInt numCells = _materialIS->size();
+ const PetscInt* cells = _materialIS->points();
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();assert(cs);
@@ -334,9 +336,9 @@
// Get cell information
PetscDM dmMesh = field->mesh().dmMesh();assert(dmMesh);
- topology::StratumIS materialIS(dmMesh, "material-id", _id);
- const PetscInt numCells = materialIS.size();
- const PetscInt* cells = materialIS.points();
+ assert(_materialIS);
+ const PetscInt numCells = _materialIS->size();
+ const PetscInt* cells = _materialIS->points();
topology::FieldBase::VectorFieldEnum fieldType = topology::FieldBase::OTHER;
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.hh 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.hh 2013-05-20 17:32:42 UTC (rev 22110)
@@ -282,6 +282,8 @@
spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer
+ topology::StratumIS* _materialIS; ///< Index set for material cells.
+
int _numPropsQuadPt; ///< Number of properties per quad point.
int _numVarsQuadPt; ///< Number of state variables per quad point.
const int _dimension; ///< Spatial dimension associated with material.
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilter.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilter.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilter.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -21,13 +21,15 @@
#include "CellFilter.hh" // Implementation of class methods
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/Stratum.hh" // USES StratumIS
#include "pylith/utils/error.h" // USES PYLITH_METHOD_BEGIN/END
// ----------------------------------------------------------------------
// Constructor
pylith::meshio::CellFilter::CellFilter(void) :
- _quadrature(0)
+ _quadrature(0),
+ _cellsIS(0)
{ // constructor
} // constructor
@@ -41,7 +43,8 @@
// ----------------------------------------------------------------------
// Copy constructor.
pylith::meshio::CellFilter::CellFilter(const CellFilter& f) :
- _quadrature(0)
+ _quadrature(0),
+ _cellsIS(0)
{ // copy constructor
PYLITH_METHOD_BEGIN;
@@ -57,6 +60,7 @@
pylith::meshio::CellFilter::deallocate(void)
{ // deallocate
delete _quadrature; _quadrature = 0;
+ delete _cellsIS; _cellsIS = 0;
} // deallocate
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilter.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilter.hh 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilter.hh 2013-05-20 17:32:42 UTC (rev 22110)
@@ -102,6 +102,8 @@
/// Quadrature associated with cells.
feassemble::Quadrature* _quadrature;
+ topology::StratumIS* _cellsIS; ///< Index set of cells.
+
}; // CellFilter
#endif // pylith_meshio_cellfilter_hh
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -24,6 +24,7 @@
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/Stratum.hh" // USES StratumIS
#include "pylith/utils/error.h" // USES PYLITH_METHOD_BEGIN/END
@@ -94,8 +95,7 @@
const scalar_array& wts = quadrature->quadWts();
PetscDM dmMesh = fieldIn.mesh().dmMesh();assert(dmMesh);
- PetscIS cellIS = NULL;
- PetscInt cStart, cEnd, numCells;
+ PetscInt cStart = 0, cEnd, numCells;
PetscErrorCode err;
if (!label) {
@@ -106,13 +106,14 @@
if (cMax >= 0) {cEnd = PetscMin(cEnd, cMax);}
numCells = cEnd - cStart;
} else {
- const PetscInt *cells = NULL;
- err = DMPlexGetStratumIS(dmMesh, label, labelId, &cellIS);PYLITH_CHECK_ERROR(err);
- err = ISGetSize(cellIS, &numCells);PYLITH_CHECK_ERROR(err);
- err = ISGetIndices(cellIS, &cells);PYLITH_CHECK_ERROR(err);
- cStart = cells[0];
- err = ISRestoreIndices(cellIS, &cells);PYLITH_CHECK_ERROR(err);
- } // if
+ if (!_cellsIS) {
+ _cellsIS = new topology::StratumIS(dmMesh, label, labelId);assert(_cellsIS);
+ } // if
+ numCells = _cellsIS->size();
+ if (numCells > 0) {
+ cStart = _cellsIS->points()[0];
+ } // if
+ } // if/else
topology::VecVisitorMesh fieldInVisitor(fieldIn);
const PetscScalar* fieldInArray = fieldInVisitor.localArray();
@@ -168,9 +169,8 @@
volume += wts[iQuad];
// Loop over cells
- if (cellIS) {
- const PetscInt *cells = NULL;
- err = ISGetIndices(cellIS, &cells);PYLITH_CHECK_ERROR(err);
+ if (_cellsIS) {
+ const PetscInt* cells = _cellsIS->points();
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt ioff = fieldInVisitor.sectionOffset(cells[c]);
assert(totalFiberDim == fieldInVisitor.sectionDof(cells[c]));
@@ -184,7 +184,6 @@
fieldAvgArray[aoff+i] += wts[iQuad] / volume * fieldInArray[ioff+iQuad*fiberDim+i];
} // for
} // for
- err = ISRestoreIndices(cellIS, &cells);PYLITH_CHECK_ERROR(err);
} else {
for(PetscInt c = cStart; c < cEnd; ++c) {
const PetscInt ioff = fieldInVisitor.sectionOffset(c);
@@ -200,7 +199,6 @@
} // for
} // for
} // if/else
- err = ISDestroy(&cellIS);PYLITH_CHECK_ERROR(err);
PetscLogFlops(numCells * numQuadPts*fiberDim*3);
PYLITH_METHOD_RETURN(*_fieldAvg);
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -286,29 +286,6 @@
if (_tstampIndex == istep)
_writeTimeStamp(t, commRank);
-#if 0 // :TODO: MATT What is this doing here?
- const int spaceDim = mesh.coordsys()->spaceDim();
- PetscInt bs;
- err = VecGetBlockSize(vector, &bs);PYLITH_CHECK_ERROR(err);
- switch (field.vectorFieldType()) {
- case pylith::topology::FieldBase::VECTOR:
- if (bs % spaceDim) {
- PYLITH_CHECK_ERROR(PETSC_ERR_ARG_WRONG);
- } // if
- break;
- case pylith::topology::FieldBase::TENSOR:
- if (bs % spaceDim) {
- PYLITH_CHECK_ERROR(PETSC_ERR_ARG_WRONG);
- } // if
- break;
- default:
- if (bs > 1) {
- PYLITH_CHECK_ERROR(PETSC_ERR_ARG_WRONG);
- } // if
- break;
- } // switch
-#endif
-
err = PetscViewerHDF5PushGroup(_viewer, "/vertex_fields");PYLITH_CHECK_ERROR(err);
err = PetscViewerHDF5SetTimestep(_viewer, istep);PYLITH_CHECK_ERROR(err);
err = VecView(vector, _viewer);PYLITH_CHECK_ERROR(err);
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2013-05-20 17:32:42 UTC (rev 22110)
@@ -25,6 +25,7 @@
#include "pylith/topology/Mesh.hh" /// USES Mesh
#include "pylith/topology/Field.hh" /// USES Field
+#include "pylith/topology/Stratum.hh" /// USES StratumIS
#include "spatialdata/geocoords/CoordSys.hh" /// USES CoordSys
@@ -83,9 +84,9 @@
// Prepare for writing files.
void
pylith::meshio::DataWriterHDF5Ext::open(const topology::Mesh& mesh,
- const int numTimeSteps,
- const char* label,
- const int labelId)
+ const int numTimeSteps,
+ const char* label,
+ const int labelId)
{ // open
PYLITH_METHOD_BEGIN;
@@ -359,47 +360,59 @@
PetscVec vector = field.vector(context);assert(vector);
err = VecView(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
- ++_datasets[field.label()].numTimeSteps;
- PetscDM dm = NULL;
- PetscSection section = NULL;
- PetscInt dof = 0, n, numLocalVertices = 0, numVertices, vStart;
- PetscIS globalVertexNumbers = NULL;
+ ExternalDataset& datasetInfo = _datasets[field.label()];
+ ++datasetInfo.numTimeSteps;
- err = VecGetDM(vector, &dm);PYLITH_CHECK_ERROR(err);
- err = DMGetDefaultSection(dm, §ion);PYLITH_CHECK_ERROR(err);assert(section);
+ // Update time stamp in "/time, if necessary.
+ if (!commRank) {
+ if (_tstampIndex+1 == datasetInfo.numTimeSteps)
+ _writeTimeStamp(t);
+ } // if
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, PETSC_NULL);PYLITH_CHECK_ERROR(err);
- err = DMPlexGetVertexNumbering(dmMesh, &globalVertexNumbers);PYLITH_CHECK_ERROR(err);
- err = ISGetSize(globalVertexNumbers, &n);PYLITH_CHECK_ERROR(err);
- if (n > 0) {
- const PetscInt *indices = NULL;
- err = ISGetIndices(globalVertexNumbers, &indices);PYLITH_CHECK_ERROR(err);
- err = PetscSectionGetDof(section, vStart, &dof);PYLITH_CHECK_ERROR(err);
- for(PetscInt v = 0; v < n; ++v) {
- if (indices[v] >= 0) ++numLocalVertices;
+ // Add dataset to HDF5 file, if necessary
+ if (createdExternalDataset) {
+ PetscDM dm = NULL;
+ PetscSection section = NULL;
+ PetscInt dof = 0, n, numLocalVertices = 0, numVertices, vStart;
+ PetscIS globalVertexNumbers = NULL;
+
+ err = VecGetDM(vector, &dm);PYLITH_CHECK_ERROR(err);
+ err = DMGetDefaultSection(dm, §ion);PYLITH_CHECK_ERROR(err);assert(section);
+
+ err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, PETSC_NULL);PYLITH_CHECK_ERROR(err);
+ err = DMPlexGetVertexNumbering(dmMesh, &globalVertexNumbers);PYLITH_CHECK_ERROR(err);
+ err = ISGetSize(globalVertexNumbers, &n);PYLITH_CHECK_ERROR(err);
+ if (n > 0) {
+ const PetscInt *indices = NULL;
+ err = ISGetIndices(globalVertexNumbers, &indices);PYLITH_CHECK_ERROR(err);
+ err = PetscSectionGetDof(section, vStart, &dof);PYLITH_CHECK_ERROR(err);
+ for(PetscInt v = 0; v < n; ++v) {
+ if (indices[v] >= 0) ++numLocalVertices;
+ } // for
+ err = ISRestoreIndices(globalVertexNumbers, &indices);PYLITH_CHECK_ERROR(err);
} // for
- err = ISRestoreIndices(globalVertexNumbers, &indices);PYLITH_CHECK_ERROR(err);
- } // for
- int fiberDimLocal = dof;
- int fiberDim = 0;
- err = MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, comm);PYLITH_CHECK_ERROR(err);
- err = MPI_Allreduce(&numLocalVertices, &numVertices, 1, MPI_INT, MPI_SUM, comm);PYLITH_CHECK_ERROR(err);
- assert(fiberDim > 0);assert(numVertices > 0);
+ int fiberDimLocal = dof;
+ int fiberDim = 0;
+ err = MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, comm);PYLITH_CHECK_ERROR(err);
+ err = MPI_Allreduce(&numLocalVertices, &numVertices, 1, MPI_INT, MPI_SUM, comm);PYLITH_CHECK_ERROR(err);
+ assert(fiberDim > 0);assert(numVertices > 0);
- if (!commRank) {
- if (createdExternalDataset) {
+ datasetInfo.numPoints = numVertices;
+ datasetInfo.fiberDim = fiberDim;
+
+ if (!commRank) {
// Add new external dataset to HDF5 file.
const int numTimeSteps = DataWriter::_numTimeSteps;
const hsize_t ndims = (numTimeSteps > 0) ? 3 : 2;
hsize_t maxDims[3];
if (3 == ndims) {
maxDims[0] = H5S_UNLIMITED;
- maxDims[1] = numVertices;
- maxDims[2] = fiberDim;
+ maxDims[1] = datasetInfo.numPoints;
+ maxDims[2] = datasetInfo.fiberDim;
} else {
- maxDims[0] = numVertices;
- maxDims[1] = fiberDim;
+ maxDims[0] = datasetInfo.numPoints;
+ maxDims[1] = datasetInfo.fiberDim;
} // else
// Create 'vertex_fields' group if necessary.
if (!_h5->hasGroup("/vertex_fields"))
@@ -408,24 +421,21 @@
_h5->createDatasetRawExternal("/vertex_fields", field.label(), _datasetFilename(field.label()).c_str(), maxDims, ndims, scalartype);
std::string fullName = std::string("/vertex_fields/") + field.label();
_h5->writeAttribute(fullName.c_str(), "vector_field_type", topology::FieldBase::vectorFieldString(field.vectorFieldType()));
- } else {
+ } // if
+ } else {
+ if (!commRank) {
// Update number of time steps in external dataset info in HDF5 file.
- const int totalNumTimeSteps =
- DataWriter::_numTimeSteps;
+ const int totalNumTimeSteps = DataWriter::_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] = numVertices;
- dims[2] = fiberDim;
+ dims[0] = datasetInfo.numTimeSteps; // update to current value
+ dims[1] = datasetInfo.numPoints;
+ dims[2] = datasetInfo.fiberDim;
_h5->extendDatasetRawExternal("/vertex_fields", field.label(), dims, ndims);
- } // if/else
- // Update time stamp in "/time, if necessary.
- if (_tstampIndex+1 == _datasets[field.label()].numTimeSteps)
- _writeTimeStamp(t);
- } // if
+ } // if
+ } // if/else
} catch (const std::exception& err) {
std::ostringstream msg;
@@ -489,64 +499,68 @@
PetscVec vector = field.vector(context);assert(vector);
err = VecView(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
- ++_datasets[field.label()].numTimeSteps;
- PetscSection section = field.petscSection();assert(section);
- PetscInt dof = 0, n, numLocalCells = 0, numCells, cellHeight, cStart, cEnd;
- PetscIS globalCellNumbers;
+ ExternalDataset& datasetInfo = _datasets[field.label()];
+ ++datasetInfo.numTimeSteps;
- err = DMPlexGetVTKCellHeight(dmMesh, &cellHeight);PYLITH_CHECK_ERROR(err);
- err = DMPlexGetHeightStratum(dmMesh, cellHeight, &cStart, &cEnd);PYLITH_CHECK_ERROR(err);
- if (label) {
- PetscIS pointIS;
+ // Update time stamp in "/time, if necessary.
+ if (!commRank) {
+ if (_tstampIndex+1 == datasetInfo.numTimeSteps)
+ _writeTimeStamp(t);
+ } // if
- DMPlexGetStratumIS(dmMesh, label, labelId, &pointIS);PYLITH_CHECK_ERROR(err);
- err = ISGetLocalSize(pointIS, &n);PYLITH_CHECK_ERROR(err);
- if (n > 0) {
- const PetscInt *indices;
- err = ISGetIndices(pointIS, &indices);PYLITH_CHECK_ERROR(err);
- for(PetscInt c = 0; c < n; ++c) {
- if ((indices[c] >= cStart) && (indices[c] < cEnd)) {
- err = PetscSectionGetDof(section, indices[0], &dof);PYLITH_CHECK_ERROR(err);
- ++numLocalCells;
- } // if
- } // for
- err = ISRestoreIndices(pointIS, &indices);PYLITH_CHECK_ERROR(err);
- } // if
- err = ISDestroy(&pointIS);PYLITH_CHECK_ERROR(err);
- } else {
- err = DMPlexGetCellNumbering(dmMesh, &globalCellNumbers);PYLITH_CHECK_ERROR(err);
- err = ISGetLocalSize(globalCellNumbers, &n);PYLITH_CHECK_ERROR(err);
- if (n > 0) {
- const PetscInt *indices = NULL;
- err = ISGetIndices(globalCellNumbers, &indices);PYLITH_CHECK_ERROR(err);
- err = PetscSectionGetDof(section, cStart, &dof);PYLITH_CHECK_ERROR(err);
- for(PetscInt v = 0; v < n; ++v) {
- if (indices[v] >= 0) ++numLocalCells;
- } // for
- err = ISRestoreIndices(globalCellNumbers, &indices);PYLITH_CHECK_ERROR(err);
- } // if
- } // if/else
- int fiberDimLocal = dof;
- int fiberDim = 0;
- MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, comm);
- err = MPI_Allreduce(&numLocalCells, &numCells, 1, MPI_INT, MPI_SUM, comm);PYLITH_CHECK_ERROR(err);
- assert(fiberDim > 0);assert(numCells > 0);
+ // Add dataset to HDF5 file, if necessary
+ if (createdExternalDataset) {
+ // Get cell information
+ PetscSection section = field.petscSection();assert(section);
+ PetscInt dof = 0, n, numLocalCells = 0, numCells, cellHeight, cStart, cEnd;
+ PetscIS globalCellNumbers;
+
+ err = DMPlexGetVTKCellHeight(dmMesh, &cellHeight);PYLITH_CHECK_ERROR(err);
+ err = DMPlexGetHeightStratum(dmMesh, cellHeight, &cStart, &cEnd);PYLITH_CHECK_ERROR(err);
+ if (label) {
+ topology::StratumIS cellsIS(dmMesh, label, labelId);
+ numCells = cellsIS.size();
+ const PetscInt* cells = (numCells > 0) ? cellsIS.points() : 0;
+ for(PetscInt c = 0; c < numCells; ++c) {
+ err = PetscSectionGetDof(section, cells[0], &dof);PYLITH_CHECK_ERROR(err);
+ if ((cells[c] >= cStart) && (cells[c] < cEnd)) {
+ ++numLocalCells;
+ } // if
+ } // for
+ } else {
+ err = DMPlexGetCellNumbering(dmMesh, &globalCellNumbers);PYLITH_CHECK_ERROR(err);
+ err = ISGetLocalSize(globalCellNumbers, &numCells);PYLITH_CHECK_ERROR(err);
+ if (numCells > 0) {
+ const PetscInt *indices = NULL;
+ err = ISGetIndices(globalCellNumbers, &indices);PYLITH_CHECK_ERROR(err);
+ err = PetscSectionGetDof(section, cStart, &dof);PYLITH_CHECK_ERROR(err);
+ for(PetscInt v = 0; v < numCells; ++v) {
+ if (indices[v] >= 0) ++numLocalCells;
+ } // for
+ err = ISRestoreIndices(globalCellNumbers, &indices);PYLITH_CHECK_ERROR(err);
+ } // if
+ } // if/else
+ int fiberDimLocal = dof;
+ int fiberDim = 0;
+ MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, comm);
+ err = MPI_Allreduce(&numLocalCells, &numCells, 1, MPI_INT, MPI_SUM, comm);PYLITH_CHECK_ERROR(err);
+ assert(fiberDim > 0);assert(numCells > 0);
- if (!commRank) {
- if (createdExternalDataset) {
- // Add new external dataset to HDF5 file.
+ datasetInfo.numPoints = numCells;
+ datasetInfo.fiberDim = fiberDim;
- const int numTimeSteps =
- DataWriter::_numTimeSteps;
+ if (!commRank) {
+ // Add new external dataset to HDF5 file.
+ const int numTimeSteps = DataWriter::_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;
+ maxDims[1] = datasetInfo.numPoints;
+ maxDims[2] = datasetInfo.fiberDim;
} else {
- maxDims[0] = numCells;
+ maxDims[0] = datasetInfo.numPoints;
maxDims[1] = fiberDim;
} // else
// Create 'cell_fields' group if necessary.
@@ -557,25 +571,20 @@
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::_numTimeSteps;
- assert(totalNumTimeSteps > 0);
- const int numTimeSteps = _datasets[field.label()].numTimeSteps;
+ } // if
+
+ } else {
+ // Update number of time steps in external dataset info in HDF5 file.
+ const int totalNumTimeSteps = DataWriter::_numTimeSteps;assert(totalNumTimeSteps > 0);
- 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);
- } // if
-
+ const hsize_t ndims = 3;
+ hsize_t dims[3];
+ dims[0] = datasetInfo.numTimeSteps; // update to current value
+ dims[1] = datasetInfo.numPoints;
+ dims[2] = datasetInfo.fiberDim;
+ _h5->extendDatasetRawExternal("/cell_fields", field.label(), dims, ndims);
+ } // if/else
+
} catch (const std::exception& err) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.hh 2013-05-19 20:48:36 UTC (rev 22109)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.hh 2013-05-20 17:32:42 UTC (rev 22110)
@@ -151,7 +151,9 @@
struct ExternalDataset {
PetscViewer viewer;
- int numTimeSteps;
+ PetscInt numTimeSteps;
+ PetscInt numPoints;
+ PetscInt fiberDim;
};
typedef std::map<std::string, ExternalDataset> dataset_type;
More information about the CIG-COMMITS
mailing list