[cig-commits] r20640 - in short/3D/PyLith/trunk: libsrc/pylith/faults libsrc/pylith/meshio libsrc/pylith/topology modulesrc/topology pylith/problems
knepley at geodynamics.org
knepley at geodynamics.org
Tue Aug 28 07:37:11 PDT 2012
Author: knepley
Date: 2012-08-28 07:37:11 -0700 (Tue, 28 Aug 2012)
New Revision: 20640
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
short/3D/PyLith/trunk/modulesrc/topology/Field.i
short/3D/PyLith/trunk/pylith/problems/Formulation.py
Log:
HDF5 now works for DMComplex, Added Python bindings for new Field interface, Fixed label propagation in CohesiveTopology
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2012-08-28 02:44:27 UTC (rev 20639)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2012-08-28 14:37:11 UTC (rev 20640)
@@ -298,12 +298,17 @@
err = DMComplexSetVTKBounds(newMesh, PETSC_DETERMINE, firstLagrangeVertexDM);CHECK_PETSC_ERROR(err);
// Renumber labels
- for(std::set<std::string>::const_iterator name = groupNames->begin(); name != groupNames->end(); ++name) {
+ std::set<std::string> names(groupNames->begin(), groupNames->end());
+ names.insert(names.begin(), "material-id");
+ for(std::set<std::string>::const_iterator name = names.begin(); name != names.end(); ++name) {
const char *lname = (*name).c_str();
IS idIS;
PetscInt n;
+ PetscBool hasLabel;
const PetscInt *ids;
+ err = DMComplexHasLabel(complexMesh, lname, &hasLabel);CHECK_PETSC_ERROR(err);
+ if (!hasLabel) continue;
err = DMComplexGetLabelSize(complexMesh, lname, &n);CHECK_PETSC_ERROR(err);
err = DMComplexGetLabelIdIS(complexMesh, lname, &idIS);CHECK_PETSC_ERROR(err);
err = ISGetIndices(idIS, &ids);CHECK_PETSC_ERROR(err);
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2012-08-28 02:44:27 UTC (rev 20639)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2012-08-28 14:37:11 UTC (rev 20640)
@@ -92,6 +92,7 @@
const std::string& filename = _hdf5Filename();
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ DM dmMesh = mesh.dmMesh();
assert(!sieveMesh.isNull());
_timesteps.clear();
@@ -107,7 +108,54 @@
err = PetscViewerHDF5Open(mesh.comm(), filename.c_str(), FILE_MODE_WRITE,
&_viewer);
CHECK_PETSC_ERROR(err);
+ ALE::Obj<numbering_type> vNumbering =
+ sieveMesh->hasLabel("censored depth") ?
+ sieveMesh->getFactory()->getNumbering(sieveMesh, "censored depth", 0) :
+ sieveMesh->getFactory()->getNumbering(sieveMesh, 0);
+ assert(!vNumbering.isNull());
+ //vNumbering->view("VERTEX NUMBERING");
+ if (dmMesh) {
+ PetscSection coordSection;
+ Vec coordinates;
+ PetscInt vStart, vEnd, vMax, verticesSize, dim, dimLocal = 0;
+
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+ assert(cs);
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordinates);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetVTKBounds(dmMesh, PETSC_NULL, &vMax);CHECK_PETSC_ERROR(err);
+ if (vMax >= 0) {vEnd = PetscMin(vEnd, vMax);}
+ for(PetscInt vertex = vStart; vertex < vEnd; ++vertex) {
+ err = PetscSectionGetDof(coordSection, vertex, &dimLocal);CHECK_PETSC_ERROR(err);
+ if (dimLocal) break;
+ }
+ err = MPI_Allreduce(&dimLocal, &dim, 1, MPIU_INT, MPI_MAX, sieveMesh->comm());CHECK_PETSC_ERROR(err);
+ verticesSize = vEnd - vStart;
+
+ PetscVec coordVec;
+ PetscScalar *coords, *c;
+
+ err = VecCreate(sieveMesh->comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, verticesSize*dim, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetBlockSize(coordVec, dim);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordinates, &c);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = 0; v < vEnd - vStart; ++v) {
+ for(PetscInt d = 0; d < dim; ++d) {
+ coords[v*dim+d] = c[v*dim+d];
+ }
+ }
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(coordinates, &c);CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) coordVec, "vertices");CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PushGroup(_viewer, "/geometry");CHECK_PETSC_ERROR(err);
+ err = VecView(coordVec, _viewer);CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
+ err = VecDestroy(&coordVec);CHECK_PETSC_ERROR(err);
+ } else {
const ALE::Obj<typename mesh_type::RealSection>& coordinatesSection =
sieveMesh->hasRealSection("coordinates_dimensioned") ?
sieveMesh->getRealSection("coordinates_dimensioned") :
@@ -123,12 +171,6 @@
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
topology::Field<mesh_type> coordinates(mesh, coordinatesSection, metadata);
coordinates.label("vertices");
- ALE::Obj<numbering_type> vNumbering =
- sieveMesh->hasLabel("censored depth") ?
- sieveMesh->getFactory()->getNumbering(sieveMesh, "censored depth", 0) :
- sieveMesh->getFactory()->getNumbering(sieveMesh, 0);
- assert(!vNumbering.isNull());
- //vNumbering->view("VERTEX NUMBERING");
coordinates.createScatterWithBC(mesh, vNumbering, context);
coordinates.scatterSectionToVector(context);
@@ -137,7 +179,62 @@
err = PetscViewerHDF5PushGroup(_viewer, "/geometry");CHECK_PETSC_ERROR(err);
err = VecView(coordinatesVector, _viewer);CHECK_PETSC_ERROR(err);
err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
+ }
+ if (dmMesh) {
+ PetscInt *cones;
+ PetscSection coneSection;
+ PetscInt vStart, vEnd, cStart, cEnd, cMax, dof, conesSize, numCorners, numCornersLocal = 0;
+
+ err = DMComplexGetConeSection(dmMesh, &coneSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCones(dmMesh, &cones);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ 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);}
+ for(PetscInt cell = cStart; cell < cEnd; ++cell) {
+ err = DMComplexGetConeSize(dmMesh, cell, &numCornersLocal);CHECK_PETSC_ERROR(err);
+ if (numCornersLocal) break;
+ }
+ err = MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, sieveMesh->comm());CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(coneSection, cEnd, &conesSize);CHECK_PETSC_ERROR(err);
+ if (label) {
+ conesSize = 0;
+ for(PetscInt cell = cStart; cell < cEnd; ++cell) {
+ PetscInt value;
+
+ err = DMComplexGetLabelValue(dmMesh, label, cell, &value);CHECK_PETSC_ERROR(err);
+ if (value == labelId) ++conesSize;
+ }
+ conesSize *= numCorners;
+ }
+
+ PetscVec cellVec;
+ PetscScalar *vertices;
+
+ err = VecCreate(sieveMesh->comm(), &cellVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(cellVec, conesSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetBlockSize(cellVec, numCorners);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(cellVec);CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) cellVec, "cells");CHECK_PETSC_ERROR(err);
+ err = VecGetArray(cellVec, &vertices);CHECK_PETSC_ERROR(err);
+ for(PetscInt cell = cStart, v = 0; cell < cEnd; ++cell) {
+ if (label) {
+ PetscInt value;
+
+ err = DMComplexGetLabelValue(dmMesh, label, cell, &value);CHECK_PETSC_ERROR(err);
+ if (value != labelId) continue;
+ }
+ for(PetscInt corner = 0; corner < numCorners; ++corner, ++v) {
+ vertices[v] = cones[cell*numCorners+corner] - vStart;
+ }
+ }
+ err = VecRestoreArray(cellVec, &vertices);CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PushGroup(_viewer, "/topology");CHECK_PETSC_ERROR(err);
+ err = VecView(cellVec, _viewer);CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PopGroup(_viewer);CHECK_PETSC_ERROR(err);
+ err = VecDestroy(&cellVec);CHECK_PETSC_ERROR(err);
+ } else {
// Account for censored cells
int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
int cellDepth = 0;
@@ -206,6 +303,7 @@
err = PetscViewerHDF5PopGroup(_viewer);CHECK_PETSC_ERROR(err);
err = VecDestroy(&elemVec);CHECK_PETSC_ERROR(err);
delete[] tmpVertices; tmpVertices = 0;
+ }
hid_t h5 = -1;
err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-08-28 02:44:27 UTC (rev 20639)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-08-28 14:37:11 UTC (rev 20640)
@@ -79,8 +79,7 @@
// Constructor with field and subfields
template<typename mesh_type, typename section_type>
pylith::topology::Field<mesh_type, section_type>::Field(const Field& src,
- int numFields,
- const int fields[]) :
+ const int fields[], int numFields) :
_mesh(src._mesh),
_section(PETSC_NULL)
{ // constructor
@@ -1432,7 +1431,7 @@
// Experimental
template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type, section_type>::addField(const std::string& name, int numComponents)
+pylith::topology::Field<mesh_type, section_type>::addField(const char *name, int numComponents)
{
// Keep track of name/components until setup
_tmpFields[name] = numComponents;
@@ -1458,7 +1457,7 @@
template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type, section_type>::updateDof(const std::string& name, const DomainEnum domain, int fiberDim)
+pylith::topology::Field<mesh_type, section_type>::updateDof(const char *name, const DomainEnum domain, int fiberDim)
{
PetscSection section;
PetscInt pStart, pEnd, f = 0;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-08-28 02:44:27 UTC (rev 20639)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-08-28 14:37:11 UTC (rev 20640)
@@ -86,7 +86,7 @@
const Metadata& metadata);
- Field(const Field& src, int numFields, const int fields[]);
+ Field(const Field& src, const int fields[], int numFields);
/// Destructor.
~Field(void);
@@ -253,11 +253,11 @@
*/
void cloneSection(const Field& src);
- void addField(const std::string& name, int numComponents);
+ void addField(const char *name, int numComponents);
void setupFields();
- void updateDof(const std::string& name, const DomainEnum domain, const int fiberDim);
+ void updateDof(const char *name, const DomainEnum domain, const int fiberDim);
/// Clear variables associated with section.
void clear(void);
Modified: short/3D/PyLith/trunk/modulesrc/topology/Field.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/Field.i 2012-08-28 02:44:27 UTC (rev 20639)
+++ short/3D/PyLith/trunk/modulesrc/topology/Field.i 2012-08-28 14:37:11 UTC (rev 20640)
@@ -161,6 +161,12 @@
*/
void cloneSection(const Field& src);
+ void addField(const char *name, int numComponents);
+
+ void setupFields();
+
+ void updateDof(const char *name, const pylith::topology::FieldBase::DomainEnum domain, const int fiberDim);
+
/// Clear variables associated with section.
void clear(void);
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2012-08-28 02:44:27 UTC (rev 20639)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2012-08-28 14:37:11 UTC (rev 20640)
@@ -535,6 +535,11 @@
if 1:
solution = self.fields.get("dispIncr(t->t+dt)")
solution.newSection(solution.VERTICES_FIELD, dimension)
+
+ solution.addField("displacement", dimension)
+ solution.setupFields()
+ solution.updateDof("displacement", solution.VERTICES_FIELD, dimension)
+
solution.vectorFieldType(solution.VECTOR)
solution.scale(lengthScale.value)
if self.splitFields():
@@ -553,9 +558,9 @@
solution.addField("pressure", 1)
solution.addField("temperature", 1)
solution.setupFields()
- solution.updateDof("displacement", solution.VERTICES_FIELD, dimension)
- solution.updateDof("pressure", solution.CELLS_FIELD, 1)
- solution.updateDof("temperature", solution.VERTICES_FIELD, 1)
+ solution.updateDof("displacement", solution.VERTICES_FIELD, dimension)
+ solution.updateDof("pressure", solution.CELLS_FIELD, 1)
+ solution.updateDof("temperature", solution.VERTICES_FIELD, 1)
solution.vectorFieldType("displacement", solution.VECTOR)
solution.scale("displacement", lengthScale.value)
for constraint in self.constraints:
More information about the CIG-COMMITS
mailing list