[cig-commits] r20708 - in short/3D/PyLith/trunk: libsrc/pylith/materials libsrc/pylith/meshio libsrc/pylith/topology pylith/problems
knepley at geodynamics.org
knepley at geodynamics.org
Mon Sep 10 12:33:25 PDT 2012
Author: knepley
Date: 2012-09-10 12:33:25 -0700 (Mon, 10 Sep 2012)
New Revision: 20708
Modified:
short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
short/3D/PyLith/trunk/pylith/problems/Formulation.py
Log:
Fixed bug in Field.updateDof(), fixed another newSection(), converting more DataWriters, CellFilterAvg, and Material
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2012-09-09 02:51:34 UTC (rev 20707)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2012-09-10 19:33:25 UTC (rev 20708)
@@ -405,16 +405,14 @@
// Allocate buffer for property field if necessary.
PetscSection fieldSection = field->petscSection();
- Vec fieldVec = field->localVector();
- bool useCurrentField = fieldSection != PETSC_NULL;
- if (fieldSection) {
- PetscInt pStart, pEnd;
+ bool useCurrentField = PETSC_FALSE;
+ PetscInt pStart, pEnd;
- err = PetscSectionGetChart(fieldSection, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- if (pEnd < 0) {
- err = DMComplexGetHeightStratum(dmMesh, 0, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetChart(fieldSection, pStart, pEnd);CHECK_PETSC_ERROR(err);
- }
+ err = PetscSectionGetChart(fieldSection, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ if (pEnd < 0) {
+ err = DMComplexGetHeightStratum(dmMesh, 0, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(fieldSection, pStart, pEnd);CHECK_PETSC_ERROR(err);
+ } else {
// check fiber dimension
PetscInt totalFiberDimCurrentLocal = 0;
PetscInt totalFiberDimCurrent = 0;
@@ -444,6 +442,7 @@
scalar_array propertiesCell(numQuadPts*numPropsQuadPt);
// Loop over cells
+ Vec fieldVec = field->localVector();
PetscScalar *fieldArray, *propertiesArray;
err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc 2012-09-09 02:51:34 UTC (rev 20707)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc 2012-09-10 19:33:25 UTC (rev 20708)
@@ -97,28 +97,33 @@
const int numQuadPts = quadrature->numQuadPts();
const scalar_array& wts = quadrature->quadWts();
- const ALE::Obj<SieveMesh>& sieveMesh = fieldIn.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
- const int depth = (0 == label) ? cellDepth : labelId;
- const std::string labelName = (0 == label) ?
- ((sieveMesh->hasLabel("censored depth")) ?
- "censored depth" : "depth") : label;
+ DM dmMesh = fieldIn.mesh().dmMesh();
+ IS cellIS = PETSC_NULL;
+ PetscInt cStart, cEnd, numCells;
+ PetscErrorCode err;
- const ALE::Obj<label_sequence>& cells =
- sieveMesh->getLabelStratum(labelName, depth);
- assert(!cells.isNull());
- const typename label_sequence::iterator cellsBegin = cells->begin();
- const typename label_sequence::iterator cellsEnd = cells->end();
+ assert(dmMesh);
+ if (!label) {
+ PetscInt cMax;
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetVTKBounds(dmMesh, &cMax, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ if (cMax >= 0) {cEnd = PetscMin(cEnd, cMax);}
+ numCells = cEnd - cStart;
+ } else {
+ err = DMComplexGetStratumIS(dmMesh, label, 1, &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ }
+
// Only processors with cells for output get the correct fiber dimension.
- const ALE::Obj<RealSection>& sectionIn = fieldIn.section();
- assert(!sectionIn.isNull());
- const int totalFiberDim = (cellsBegin != cellsEnd) ?
- sectionIn->getFiberDimension(*cellsBegin) : 0;
+ PetscSection sectionIn = fieldIn.petscSection();
+ Vec vecIn = fieldIn.localVector();
+ PetscInt totalFiberDim = 0;
+
+ assert(sectionIn);
+ if (numCells) {err = PetscSectionGetDof(sectionIn, cStart, &totalFiberDim);CHECK_PETSC_ERROR(err);}
const int fiberDim = totalFiberDim / numQuadPts;
assert(fiberDim * numQuadPts == totalFiberDim);
-
// Allocate field if necessary
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("OutputFields");
@@ -127,7 +132,7 @@
assert(0 != _fieldAvg);
_fieldAvg->newSection(fieldIn, fiberDim);
_fieldAvg->allocate();
- } else if (_fieldAvg->sectionSize() != cells->size()*fiberDim) {
+ } else if (_fieldAvg->sectionSize() != numCells*fiberDim) {
_fieldAvg->newSection(fieldIn, fiberDim);
_fieldAvg->allocate();
} // else
@@ -158,32 +163,63 @@
assert(0);
throw std::logic_error("Bad vector field type for CellFilterAvg.");
} // switch
- const ALE::Obj<RealSection>& sectionAvg = _fieldAvg->section();
+ PetscSection sectionAvg = _fieldAvg->petscSection();
+ Vec vecAvg = _fieldAvg->localVector();
_fieldAvg->label(fieldIn.label());
_fieldAvg->scale(fieldIn.scale());
_fieldAvg->addDimensionOkay(true);
- scalar_array fieldAvgCell(fiberDim);
- PylithScalar scalar = 0.0;
+ PylithScalar volume = 0.0;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- scalar += wts[iQuad];
+ volume += wts[iQuad];
// Loop over cells
- for (typename label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- const PylithScalar* values = sectionIn->restrictPoint(*c_iter);
- assert(totalFiberDim == sectionIn->getFiberDimension(*c_iter));
+ PetscScalar *arrayIn, *arrayAvg;
+
+ err = VecGetArray(vecIn, &arrayIn);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(vecAvg, &arrayAvg);CHECK_PETSC_ERROR(err);
+ if (cellIS) {
+ const PetscInt *cells;
+
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = 0; c < numCells; ++c) {
+ PetscInt dof, off, adof, aoff;
- fieldAvgCell = 0.0;
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- for (int i=0; i < fiberDim; ++i)
- fieldAvgCell[i] += wts[iQuad] / scalar * values[iQuad*fiberDim+i];
+ err = PetscSectionGetDof(sectionIn, cells[c], &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(sectionIn, cells[c], &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(sectionAvg, cells[c], &adof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(sectionAvg, cells[c], &aoff);CHECK_PETSC_ERROR(err);
+ assert(totalFiberDim == dof);
+ assert(fiberDim == adof);
+ for(int i = 0; i < adof; ++i) {
+ arrayAvg[aoff+i] = 0.0;
+ for(int iQuad = 0; iQuad < numQuadPts; ++iQuad)
+ arrayAvg[aoff+i] += wts[iQuad] / volume * arrayIn[off+iQuad*fiberDim+i];
+ }
+ }
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
+ } else {
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ PetscInt dof, off, adof, aoff;
+
+ err = PetscSectionGetDof(sectionIn, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(sectionIn, c, &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(sectionAvg, c, &adof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(sectionAvg, c, &aoff);CHECK_PETSC_ERROR(err);
+ assert(totalFiberDim == dof);
+ assert(fiberDim == adof);
+ for(int i = 0; i < adof; ++i) {
+ arrayAvg[aoff+i] = 0.0;
+ for(int iQuad = 0; iQuad < numQuadPts; ++iQuad)
+ arrayAvg[aoff+i] += wts[iQuad] / volume * arrayIn[off+iQuad*fiberDim+i];
+ }
+ } // for
+ }
+ err = VecRestoreArray(vecIn, &arrayIn);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(vecAvg, &arrayAvg);CHECK_PETSC_ERROR(err);
+ PetscLogFlops(numCells * numQuadPts*fiberDim*3);
- sectionAvg->updatePoint(*c_iter, &fieldAvgCell[0]);
- } // for
- PetscLogFlops( cells->size() * numQuadPts*fiberDim*3 );
-
return *_fieldAvg;
} // filter
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2012-09-09 02:51:34 UTC (rev 20707)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2012-09-10 19:33:25 UTC (rev 20708)
@@ -93,9 +93,11 @@
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ assert(dmMesh);
+ PetscMPIInt commRank;
+ PetscErrorCode err = MPI_Comm_rank(((PetscObject) dmMesh)->comm, &commRank);CHECK_PETSC_ERROR(err);
- const int commRank = sieveMesh->commRank();
if (!commRank) {
_h5->open(_hdf5Filename().c_str(), H5F_ACC_TRUNC);
@@ -106,7 +108,6 @@
_tstampIndex = 0;
PetscViewer binaryViewer;
- PetscErrorCode err = 0;
const hid_t scalartype = (sizeof(double) == sizeof(PylithScalar)) ?
H5T_IEEE_F64BE : H5T_IEEE_F32BE;
@@ -137,7 +138,7 @@
assert(coordinatesVector);
const std::string& filenameVertices = _datasetFilename("vertices");
- err = PetscViewerBinaryOpen(sieveMesh->comm(), filenameVertices.c_str(),
+ err = PetscViewerBinaryOpen(((PetscObject) dmMesh)->comm, filenameVertices.c_str(),
FILE_MODE_WRITE,
&binaryViewer);
CHECK_PETSC_ERROR(err);
@@ -163,7 +164,7 @@
int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
int cellDepth = 0;
err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
- sieveMesh->comm());CHECK_PETSC_ERROR(err);
+ ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
((sieveMesh->hasLabel("censored depth")) ?
@@ -180,7 +181,7 @@
numCornersLocal = sieveMesh->getNumCellCorners(*cells->begin());
int numCorners = numCornersLocal;
err = MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPI_INT, MPI_MAX,
- sieveMesh->comm()); CHECK_PETSC_ERROR(err);
+ ((PetscObject) dmMesh)->comm); CHECK_PETSC_ERROR(err);
PylithScalar* tmpVertices = 0;
const int ncells = cNumbering->getLocalSize();
@@ -220,13 +221,13 @@
} // if
PetscVec elemVec;
- err = VecCreateMPIWithArray(sieveMesh->comm(), numCorners, conesSize, PETSC_DETERMINE,
+ err = VecCreateMPIWithArray(((PetscObject) dmMesh)->comm, numCorners, conesSize, PETSC_DETERMINE,
tmpVertices, &elemVec); CHECK_PETSC_ERROR(err);
err = PetscObjectSetName((PetscObject) elemVec,
"cells");CHECK_PETSC_ERROR(err);
const std::string& filenameCells = _datasetFilename("cells");
- err = PetscViewerBinaryOpen(sieveMesh->comm(), filenameCells.c_str(),
+ err = PetscViewerBinaryOpen(((PetscObject) dmMesh)->comm, filenameCells.c_str(),
FILE_MODE_WRITE,
&binaryViewer);
CHECK_PETSC_ERROR(err);
@@ -302,11 +303,13 @@
try {
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
- const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = field.mesh().sieveMesh();
assert(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ assert(dmMesh);
+ PetscMPIInt commRank;
+ PetscErrorCode err = MPI_Comm_rank(((PetscObject) dmMesh)->comm, &commRank);CHECK_PETSC_ERROR(err);
- const int commRank = sieveMesh->commRank();
-
const std::string labelName =
(sieveMesh->hasLabel("censored depth")) ? "censored depth" : "depth";
ALE::Obj<numbering_type> vNumbering =
@@ -320,7 +323,6 @@
assert(vector);
PetscViewer binaryViewer;
- PetscErrorCode err = 0;
const hid_t scalartype = (sizeof(double) == sizeof(PylithScalar)) ?
H5T_IEEE_F64BE : H5T_IEEE_F32BE;
@@ -330,7 +332,7 @@
if (_datasets.find(field.label()) != _datasets.end()) {
binaryViewer = _datasets[field.label()].viewer;
} else {
- err = PetscViewerBinaryOpen(sieveMesh->comm(),
+ err = PetscViewerBinaryOpen(((PetscObject) dmMesh)->comm,
_datasetFilename(field.label()).c_str(),
FILE_MODE_WRITE, &binaryViewer);
CHECK_PETSC_ERROR(err);
@@ -349,63 +351,60 @@
CHECK_PETSC_ERROR(err);
++_datasets[field.label()].numTimeSteps;
- const ALE::Obj<typename mesh_type::RealSection>& section =
- field.section();
- assert(!section.isNull());
+ PetscSection section = field.petscSection();
+ PetscInt dof = 0;
+ assert(section);
assert(!sieveMesh->getLabelStratum(labelName, 0).isNull());
- int fiberDimLocal =
- (sieveMesh->getLabelStratum(labelName, 0)->size() > 0) ?
- section->getFiberDimension(*sieveMesh->getLabelStratum(labelName,
- 0)->begin()) : 0;
+ if (sieveMesh->getLabelStratum(labelName, 0)->size() > 0) {
+ err = PetscSectionGetDof(section, *sieveMesh->getLabelStratum(labelName, 0)->begin(), &dof);CHECK_PETSC_ERROR(err);
+ }
+ int fiberDimLocal = dof;
int fiberDim = 0;
- MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX,
- field.mesh().comm());
+ MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, ((PetscObject) dmMesh)->comm);
assert(fiberDim > 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] = vNumbering->getGlobalSize();
- maxDims[2] = fiberDim;
- } else {
- maxDims[0] = vNumbering->getGlobalSize();
- maxDims[1] = fiberDim;
- } // else
- // Create 'vertex_fields' group if necessary.
- if (!_h5->hasGroup("/vertex_fields"))
- _h5->createGroup("/vertex_fields");
+ // 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] = vNumbering->getGlobalSize();
+ maxDims[2] = fiberDim;
+ } else {
+ maxDims[0] = vNumbering->getGlobalSize();
+ maxDims[1] = fiberDim;
+ } // else
+ // Create 'vertex_fields' group if necessary.
+ if (!_h5->hasGroup("/vertex_fields"))
+ _h5->createGroup("/vertex_fields");
- _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()));
-
+ _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 {
- // 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] = vNumbering->getGlobalSize();
- dims[2] = fiberDim;
- _h5->extendDatasetRawExternal("/vertex_fields", field.label(),
- dims, ndims);
+ const hsize_t ndims = 3;
+ hsize_t dims[3];
+ dims[0] = numTimeSteps; // update to current value
+ dims[1] = vNumbering->getGlobalSize();
+ dims[2] = 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);
+ _writeTimeStamp(t);
} // if
} catch (const std::exception& err) {
@@ -437,18 +436,18 @@
try {
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
- PetscErrorCode err = 0;
- const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
- field.mesh().sieveMesh();
+ const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = field.mesh().sieveMesh();
assert(!sieveMesh.isNull());
+ DM dmMesh = field.mesh().dmMesh();
+ assert(dmMesh);
+ PetscMPIInt commRank;
+ PetscErrorCode err = MPI_Comm_rank(((PetscObject) dmMesh)->comm, &commRank);CHECK_PETSC_ERROR(err);
- const int commRank = sieveMesh->commRank();
-
int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
int cellDepth = 0;
err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
- sieveMesh->comm());CHECK_PETSC_ERROR(err);
+ ((PetscObject) dmMesh)->comm);CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
((sieveMesh->hasLabel("censored depth")) ?
@@ -472,7 +471,7 @@
if (_datasets.find(field.label()) != _datasets.end()) {
binaryViewer = _datasets[field.label()].viewer;
} else {
- err = PetscViewerBinaryOpen(sieveMesh->comm(),
+ err = PetscViewerBinaryOpen(((PetscObject) dmMesh)->comm,
_datasetFilename(field.label()).c_str(),
FILE_MODE_WRITE, &binaryViewer);
CHECK_PETSC_ERROR(err);
@@ -491,19 +490,18 @@
CHECK_PETSC_ERROR(err);
++_datasets[field.label()].numTimeSteps;
- assert(!sieveMesh->getLabelStratum(labelName, depth).isNull());
- const ALE::Obj<typename mesh_type::RealSection>& section = field.section();
- assert(!section.isNull());
-
- int fiberDimLocal =
- (sieveMesh->getLabelStratum(labelName, depth)->size() > 0) ?
- section->getFiberDimension(*sieveMesh->getLabelStratum(labelName, depth)->begin()) : 0;
+ PetscSection section = field.petscSection();
+ PetscInt dof = 0;
+ 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);
+ }
+ int fiberDimLocal = dof;
int fiberDim = 0;
- MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX,
- field.mesh().comm());
+ MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, ((PetscObject) dmMesh)->comm);
assert(fiberDim > 0);
-
if (!commRank) {
if (createdExternalDataset) {
// Add new external dataset to HDF5 file.
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc 2012-09-09 02:51:34 UTC (rev 20707)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc 2012-09-10 19:33:25 UTC (rev 20708)
@@ -253,13 +253,15 @@
sieveMesh->getFactory()->getNumbering(sieveMesh, labelName, 0);
assert(!numbering.isNull());
- const ALE::Obj<RealSection>& section = field.section();
- assert(!section.isNull());
+ PetscSection sectionP = field.petscSection();
+ PetscInt dof = 0;
+
+ assert(sectionP);
assert(!sieveMesh->getLabelStratum(labelName, 0).isNull());
-
- int fiberDimLocal =
- (sieveMesh->getLabelStratum(labelName, 0)->size() > 0) ?
- section->getFiberDimension(*sieveMesh->getLabelStratum(labelName, 0)->begin()) : 0;
+ if (sieveMesh->getLabelStratum(labelName, 0)->size() > 0) {
+ PetscErrorCode err = PetscSectionGetDof(sectionP, *sieveMesh->getLabelStratum(labelName, 0)->begin(), &dof);CHECK_PETSC_ERROR(err);
+ }
+ int fiberDimLocal = dof;
int fiberDim = 0;
MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX,
field.mesh().comm());
@@ -275,7 +277,7 @@
_wroteVertexHeader = true;
} // if
- err = VTKViewer::writeField(section, field.label(), fiberDim, numbering,
+ err = VTKViewer::writeField(field.section(), field.label(), fiberDim, numbering,
_viewer, enforceDim, _precision);
CHECK_PETSC_ERROR(err);
}
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.cc 2012-09-09 02:51:34 UTC (rev 20707)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.cc 2012-09-10 19:33:25 UTC (rev 20708)
@@ -196,7 +196,11 @@
(0 == _cellFilter) ? field : _cellFilter->filter(field, label, labelId);
field_type& fieldDimensioned = _dimension(fieldFiltered);
- _writer->writeCellField(t, fieldDimensioned, label, labelId);
+ try {
+ _writer->writeCellField(t, fieldDimensioned, label, labelId);
+ } catch(std::runtime_error e) {
+ std::cout << "ERROR: " << e.what() << std::endl<<std::endl<<std::endl;
+ }
} // appendCellField
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-09 02:51:34 UTC (rev 20707)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-10 19:33:25 UTC (rev 20708)
@@ -377,6 +377,18 @@
++c_iter)
if (srcSection->getFiberDimension(*c_iter) > 0)
_section->setFiberDimension(*c_iter, fiberDim);
+
+ if (_dm) {
+ PetscSection s;
+ PetscInt pStart = srcSection->getChart().min(), pEnd = srcSection->getChart().max();
+ PetscErrorCode err;
+
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, pStart, pEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ err = PetscSectionSetDof(s, p, fiberDim);CHECK_PETSC_ERROR(err);
+ }
+ }
} // if
logger.stagePop();
@@ -1577,8 +1589,10 @@
}
err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
assert(section);
- for(map_type::const_iterator f_iter = _metadata.begin(); f_iter != _metadata.end(); ++f_iter, ++f) {
+ for(map_type::const_iterator f_iter = _metadata.begin(); f_iter != _metadata.end(); ++f_iter) {
if (f_iter->first == name) break;
+ if (f_iter->first == "default") continue;
+ ++f;
}
assert(f < _metadata.size());
for(PetscInt p = pStart; p < pEnd; ++p) {
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-09-09 02:51:34 UTC (rev 20707)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-09-10 19:33:25 UTC (rev 20708)
@@ -29,6 +29,8 @@
pylith::topology::Field<mesh_type, section_type>::section(void) const {
std::ostringstream msg;
+ char name[1024];
+ PetscGetProgramName(name, 1024);
msg << "Sieve Sections are no longer supported: " << const_cast<Field*>(this)->_metadata["default"].label << "\n"
<< " Destination section:\n"
<< " space dim: " << spaceDim() << "\n"
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2012-09-09 02:51:34 UTC (rev 20707)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2012-09-10 19:33:25 UTC (rev 20708)
@@ -534,10 +534,10 @@
lengthScale = normalizer.lengthScale()
if 1:
solution = self.fields.get("dispIncr(t->t+dt)")
- solution.newSection(solution.VERTICES_FIELD, dimension)
-
solution.addField("displacement", dimension)
solution.setupFields()
+
+ solution.newSection(solution.VERTICES_FIELD, dimension)
solution.updateDof("displacement", solution.VERTICES_FIELD, dimension)
solution.vectorFieldType(solution.VECTOR)
More information about the CIG-COMMITS
mailing list