[cig-commits] r21419 - in short/3D/PyLith/trunk: libsrc/pylith/meshio libsrc/pylith/topology tests_auto/1d/line2
knepley at geodynamics.org
knepley at geodynamics.org
Thu Feb 28 15:51:31 PST 2013
Author: knepley
Date: 2013-02-28 15:51:31 -0800 (Thu, 28 Feb 2013)
New Revision: 21419
Modified:
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/libsrc/pylith/topology/MeshOps.cc
short/3D/PyLith/trunk/tests_auto/1d/line2/extensiondisp_np2.cfg
Log:
Fixed bug in HDF% vertex output, fixed some MeshOps
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2013-02-28 23:43:25 UTC (rev 21418)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2013-02-28 23:51:31 UTC (rev 21419)
@@ -126,47 +126,22 @@
err = VecView(coordVec, _viewer);CHECK_PETSC_ERROR(err);
err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
#else
- PetscSection coordSection;
- Vec coordinates;
- PetscReal lengthScale;
- PetscInt vStart, vEnd, vMax, verticesSize, dim, dimLocal = 0;
-
- /* TODO Get rid of this and use the createScatterWithBC(numbering) code */
- err = DMPlexGetScale(dmMesh, PETSC_UNIT_LENGTH, &lengthScale);CHECK_PETSC_ERROR(err);
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ DM dmCoord;
+ Vec coordinates;
+ topology::FieldBase::Metadata metadata;
+ metadata.label = "vertices";
+ metadata.vectorFieldType = topology::FieldBase::VECTOR;
+ err = DMGetCoordinateDM(dmMesh, &dmCoord);CHECK_PETSC_ERROR(err);
+ err = PetscObjectReference((PetscObject) dmCoord);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(dmMesh, &coordinates);CHECK_PETSC_ERROR(err);
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
- err = DMPlexGetHybridBounds(dmMesh, PETSC_NULL, PETSC_NULL, 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, mesh.comm());CHECK_PETSC_ERROR(err);
- verticesSize = vEnd - vStart;
-
- PetscVec coordVec;
- PetscScalar *coords, *c;
-
- err = VecCreate(mesh.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 = VecScale(coordVec, lengthScale);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) coordVec, "vertices");CHECK_PETSC_ERROR(err);
+ topology::Field<mesh_type> field(mesh, dmCoord, coordinates, metadata);
+ field.createScatterWithBC(mesh, "", 0, metadata.label.c_str());
+ field.scatterSectionToVector(metadata.label.c_str());
+ PetscVec vector = field.vector(metadata.label.c_str());
+ assert(vector);
err = PetscViewerHDF5PushGroup(_viewer, "/geometry");CHECK_PETSC_ERROR(err);
- err = VecView(coordVec, _viewer);CHECK_PETSC_ERROR(err);
+ err = VecView(vector, _viewer);CHECK_PETSC_ERROR(err);
err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
- err = VecDestroy(&coordVec);CHECK_PETSC_ERROR(err);
#endif
#if 0
const ALE::Obj<typename mesh_type::RealSection>& coordinatesSection =
@@ -191,7 +166,7 @@
err = VecView(coordinatesVector, _viewer);CHECK_PETSC_ERROR(err);
err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
#endif
- PetscInt cStart, cEnd, cMax, dof, conesSize, numCorners, numCornersLocal = 0;
+ PetscInt vStart, vEnd, cStart, cEnd, cMax, dof, conesSize, numCorners, numCornersLocal = 0;
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2013-02-28 23:43:25 UTC (rev 21418)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2013-02-28 23:51:31 UTC (rev 21419)
@@ -85,6 +85,23 @@
} // constructor
// ----------------------------------------------------------------------
+// Constructor with mesh, DM, local data, and metadata
+template<typename mesh_type>
+pylith::topology::Field<mesh_type>::Field(const mesh_type& mesh, DM dm, Vec localVec, const Metadata& metadata) :
+ _mesh(mesh),
+ _dm(dm)
+{ // constructor
+ PetscErrorCode err;
+
+ _metadata["default"] = metadata;
+ err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
+ err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) _localVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+ err = VecCopy(localVec, _localVec);CHECK_PETSC_ERROR(err);
+} // constructor
+
+// ----------------------------------------------------------------------
// Constructor with field and subfields
template<typename mesh_type>
pylith::topology::Field<mesh_type>::Field(const Field& src, const int fields[], int numFields) :
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2013-02-28 23:43:25 UTC (rev 21418)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2013-02-28 23:51:31 UTC (rev 21419)
@@ -73,13 +73,12 @@
*/
Field(const mesh_type& mesh, DM dm, const Metadata& metadata);
- /** Constructor with mesh, section, and metadata.
+ /** Constructor with mesh, DM, local data, and metadata
*
* @param mesh Finite-element mesh.
*/
- Field(const mesh_type& mesh, const Metadata& metadata);
+ Field(const mesh_type& mesh, DM dm, Vec localVec, const Metadata& metadata);
-
Field(const Field& src, const int fields[], int numFields);
/// Destructor.
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.cc 2013-02-28 23:43:25 UTC (rev 21418)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOps.cc 2013-02-28 23:51:31 UTC (rev 21419)
@@ -36,15 +36,12 @@
pylith::topology::MeshOps::numMaterialCells(const Mesh& mesh,
int materialId)
{ // numMaterialCells
- int ncells = 0;
+ PetscInt ncells = 0;
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- if (!sieveMesh.isNull()) {
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- if (!cells.isNull())
- ncells = cells->size();
- } // if
+ DM dmMesh = mesh.dmMesh();
+ assert(dmMesh);
+ PetscErrorCode err = DMPlexGetStratumSize(dmMesh, "material-id", materialId, &ncells);CHECK_PETSC_ERROR(err);
+ return ncells;
} // numMaterialCells
@@ -54,11 +51,9 @@
int* const materialIds,
const int numMaterials)
{ // checkMaterialIds
- typedef Mesh::SieveMesh SieveMesh;
+ assert((!numMaterials && !materialIds) || (numMaterials && materialIds));
+ PetscErrorCode err;
- assert( (0 == numMaterials && 0 == materialIds) ||
- (0 < numMaterials && 0 != materialIds) );
-
// Create map with indices for each material
std::map<int, int> materialIndex;
for (int i=0; i < numMaterials; ++i)
@@ -67,45 +62,42 @@
int_array matCellCounts(numMaterials);
matCellCounts = 0;
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- const ALE::Obj<SieveMesh::label_type>& materialsLabel =
- sieveMesh->getLabel("material-id");
+ DM dmMesh = mesh.dmMesh();
+ PetscInt cStart, cEnd;
- //mesh.view("===== MESH BEFORE CHECKING MATERIALS =====");
- //materialsLabel->view("material-id");
+ assert(dmMesh);
+ err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ DMLabel materialsLabel;
+ err = DMPlexGetLabel(dmMesh, "material-id", &materialsLabel);CHECK_PETSC_ERROR(err);
+ assert(materialsLabel);
- int* matBegin = materialIds;
- int* matEnd = materialIds + numMaterials;
+ int *matBegin = materialIds;
+ int *matEnd = materialIds + numMaterials;
std::sort(matBegin, matEnd);
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- const int cellId = sieveMesh->getValue(materialsLabel, *c_iter);
- const int* result = std::find(matBegin, matEnd, cellId);
+ for (PetscInt c = cStart; c < cEnd; ++c) {
+ PetscInt matId;
+
+ err = DMLabelGetValue(materialsLabel, c, &matId);CHECK_PETSC_ERROR(err);
+ const int *result = std::find(matBegin, matEnd, matId);
if (result == matEnd) {
std::ostringstream msg;
- msg << "Material id '" << cellId << "' for cell '" << *c_iter
- << "' does not match the id of any available materials or interfaces.";
+ msg << "Material id '" << matId << "' for cell '" << c
+ << "' does not match the id of any available materials or interfaces.";
throw std::runtime_error(msg.str());
} // if
- const int matIndex = materialIndex[cellId];
+ const int matIndex = materialIndex[matId];
assert(0 <= matIndex && matIndex < numMaterials);
++matCellCounts[matIndex];
} // for
// Make sure each material has
int_array matCellCountsAll(matCellCounts.size());
- MPI_Allreduce(&matCellCounts[0], &matCellCountsAll[0],
- matCellCounts.size(), MPI_INT, MPI_SUM, mesh.comm());
+ err = MPI_Allreduce(&matCellCounts[0], &matCellCountsAll[0],
+ matCellCounts.size(), MPI_INT, MPI_SUM, mesh.comm());CHECK_PETSC_ERROR(err);
for (int i=0; i < numMaterials; ++i) {
- const int matId = materialIds[i];
+ const int matId = materialIds[i];
const int matIndex = materialIndex[matId];
assert(0 <= matIndex && matIndex < numMaterials);
if (matCellCountsAll[matIndex] <= 0) {
Modified: short/3D/PyLith/trunk/tests_auto/1d/line2/extensiondisp_np2.cfg
===================================================================
--- short/3D/PyLith/trunk/tests_auto/1d/line2/extensiondisp_np2.cfg 2013-02-28 23:43:25 UTC (rev 21418)
+++ short/3D/PyLith/trunk/tests_auto/1d/line2/extensiondisp_np2.cfg 2013-02-28 23:51:31 UTC (rev 21419)
@@ -36,6 +36,9 @@
#formulation.solver = pylith.problems.SolverNonlinear
+[extensiondisp_np2.timedependent.formulation]
+matrix_type = aij
+
[extensiondisp_np2.timedependent.formulation.time_step]
total_time = 0.0*s
More information about the CIG-COMMITS
mailing list