[cig-commits] r22043 - in short/3D/PyLith/trunk/libsrc/pylith: meshio topology
knepley at geodynamics.org
knepley at geodynamics.org
Fri May 10 21:08:39 PDT 2013
Author: knepley
Date: 2013-05-10 21:08:39 -0700 (Fri, 10 May 2013)
New Revision: 22043
Modified:
short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc
Log:
Removed sieve from MeshBuilder, removed sieve from nondimensionalize in Mesh
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc 2013-05-11 03:43:34 UTC (rev 22042)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc 2013-05-11 04:08:39 UTC (rev 22043)
@@ -26,16 +26,11 @@
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-#include "Selection.hh" // USES boundary()
-
#include <cassert> // USES assert()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-
-// ----------------------------------------------------------------------
// Set vertices and cells in mesh.
void
pylith::meshio::MeshBuilder::buildMesh(topology::Mesh* mesh,
@@ -57,10 +52,9 @@
PetscInt dim = meshDim;
PetscMPIInt commRank = mesh->commRank();
PetscErrorCode err;
- const bool useSieve = true; // :TODO: Remove sieve
{ // Check to make sure every vertex is in at least one cell.
- // This is required by Sieve
+ // This is required by PETSc
std::vector<bool> vertexInCell(numVertices, false);
const int size = cells.size();
for (int i=0; i < size; ++i)
@@ -77,82 +71,6 @@
} // if
} // check
- /* Sieve */
- if (useSieve) {
- MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
- mesh->createSieveMesh(dim);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
- assert(!sieveMesh.isNull());
-
- ALE::Obj<SieveMesh::sieve_type> sieve = new SieveMesh::sieve_type(mesh->comm());
- sieveMesh->setSieve(sieve);
-
- if (0 == commRank || isParallel) {
- assert(coordinates->size() == numVertices*spaceDim);
- assert(cells.size() == numCells*numCorners);
- if (!interpolate) {
- // Create the ISieve
- sieve->setChart(SieveMesh::sieve_type::chart_type(0, numCells+numVertices));
- // Set cone and support sizes
- for (int c = 0; c < numCells; ++c)
- sieve->setConeSize(c, numCorners);
- sieve->symmetrizeSizes(numCells, numCorners, const_cast<int*>(&cells[0]), numCells);
- // Allocate point storage
- sieve->allocate();
- // Fill up cones
- int *cone = new int[numCorners];
- int *coneO = new int[numCorners];
- for (int v = 0; v < numCorners; ++v)
- coneO[v] = 1;
- for (int c = 0; c < numCells; ++c) {
- for (int v = 0; v < numCorners; ++v)
- cone[v] = cells[c*numCorners+v]+numCells;
- sieve->setCone(cone, c);
- sieve->setConeOrientation(coneO, c);
- } // for
- delete[] cone; cone = 0;
- delete[] coneO; coneO = 0;
- // Symmetrize to fill up supports
- sieve->symmetrize();
- } else {
- // Same old thing
- ALE::Obj<SieveFlexMesh::sieve_type> s = new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
- ALE::Obj<SieveFlexMesh::arrow_section_type> orientation = new SieveFlexMesh::arrow_section_type(sieve->comm(), sieve->debug());
-
- s->setDebug(2);
- ALE::SieveBuilder<SieveFlexMesh>::buildTopology(s, meshDim, numCells, const_cast<int*>(&cells[0]), numVertices, interpolate,
- numCorners, -1, orientation);
- std::map<SieveFlexMesh::point_type,SieveFlexMesh::point_type> renumbering;
- ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
- } // if/else
- if (!interpolate) {
- // Optimized stratification
- const ALE::Obj<SieveMesh::label_type>& height = sieveMesh->createLabel("height");
- const ALE::Obj<SieveMesh::label_type>& depth = sieveMesh->createLabel("depth");
-
- for(int c = 0; c < numCells; ++c) {
- height->setCone(0, c);
- depth->setCone(1, c);
- } // for
- for(int v = numCells; v < numCells+numVertices; ++v) {
- height->setCone(1, v);
- depth->setCone(0, v);
- } // for
- sieveMesh->setHeight(1);
- sieveMesh->setDepth(1);
- } else {
- sieveMesh->stratify();
- } // if/else
- } else {
- sieveMesh->getSieve()->setChart(SieveMesh::sieve_type::chart_type());
- sieveMesh->getSieve()->allocate();
- sieveMesh->stratify();
- } // if/else
-
- ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim, &(*coordinates)[0]);
- sieveMesh->getFactory()->clear();
- } // if
-
/* DMPlex */
PetscDM dmMesh;
PetscBool pInterpolate = interpolate ? PETSC_TRUE : PETSC_FALSE;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc 2013-05-11 03:43:34 UTC (rev 22042)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc 2013-05-11 04:08:39 UTC (rev 22043)
@@ -35,11 +35,6 @@
#include <cassert> // USES assert()
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::IntSection IntSection;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::topology::Distributor::Distributor(void)
{ // constructor
@@ -62,40 +57,7 @@
assert(newMesh);
newMesh->coordsys(origMesh.coordsys());
-
- const int commRank = origMesh.commRank();
journal::info_t info("distributor");
-
- if (0 == strcasecmp(partitioner, "")) {
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Distributing mesh using dumb partitioner." << journal::endl;
- } // if
- _distribute<ALE::DistributionNew<SieveMesh> >(newMesh, origMesh);
-#if defined(PETSC_HAVE_CHACO)
- } else if (0 == strcasecmp(partitioner, "chaco")) {
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Distributing mesh using 'chaco' partitioner." << journal::endl;
- } // if
- _distribute<ALE::DistributionNew<SieveMesh, ALE::Partitioner<ALE::Chaco::Partitioner<> > > >(newMesh, origMesh);
-#endif
-#if defined(PETSC_HAVE_PARMETIS)
- } else if (0 == strcasecmp(partitioner, "parmetis")) {
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Distributing mesh using 'parmetis' partitioner." << journal::endl;
- } // if
- _distribute<ALE::DistributionNew<SieveMesh, ALE::Partitioner<ALE::ParMetis::Partitioner<> > > >(newMesh, origMesh);
-#endif
- } else {
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Unknown partitioner '" << partitioner
- << "', distribution mesh using dumb partitioner." << journal::endl;
- } // if
- _distribute<ALE::DistributionNew<SieveMesh> >(newMesh, origMesh);
- } // else
PetscDM newDM = NULL;
PetscErrorCode err = DMPlexDistribute(origMesh.dmMesh(), partitioner, 0, &newDM);PYLITH_CHECK_ERROR(err);
@@ -157,258 +119,4 @@
PYLITH_METHOD_END;
} // write
-// ----------------------------------------------------------------------
-template<typename DistributionType>
-void
-pylith::topology::Distributor::_distribute(topology::Mesh* const newMesh,
- const topology::Mesh& origMesh)
-{ // distribute
- PYLITH_METHOD_BEGIN;
-
- typedef typename SieveMesh::point_type point_type;
- typedef typename DistributionType::partitioner_type partitioner_type;
- typedef typename DistributionType::partition_type partition_type;
-
- journal::info_t info("distributor");
-
- const int commRank = origMesh.commRank();
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Partitioning and distributing the mesh." << journal::endl;
- } // if
-
- ALE::Obj<SieveMesh>& newSieveMesh = newMesh->sieveMesh();
- assert(!newSieveMesh.isNull());
- const ALE::Obj<SieveMesh>& origSieveMesh = origMesh.sieveMesh();
- assert(!origSieveMesh.isNull());
-
- const ALE::Obj<SieveMesh::sieve_type> newSieve =
- new SieveMesh::sieve_type(origSieveMesh->comm(), origSieveMesh->debug());
- assert(!newSieve.isNull());
- const ALE::Obj<SieveMesh::send_overlap_type> sendMeshOverlap =
- new SieveMesh::send_overlap_type(origSieveMesh->comm(),
- origSieveMesh->debug());
- assert(!sendMeshOverlap.isNull());
- const ALE::Obj<SieveMesh::recv_overlap_type> recvMeshOverlap =
- new SieveMesh::recv_overlap_type(origSieveMesh->comm(),
- origSieveMesh->debug());
- assert(!recvMeshOverlap.isNull());
-
- newSieveMesh = new SieveMesh(origSieveMesh->comm(),
- origSieveMesh->getDimension(),
- origSieveMesh->debug());
- assert(!newSieveMesh.isNull());
- newSieveMesh->setSieve(newSieve);
- // IMESH_TODO This might be unnecessary, since the overlap for
- // submeshes is just the restriction of the overlaps
- // std::map<point_type,point_type> renumbering;
- SieveMesh::renumbering_type& renumbering = newSieveMesh->getRenumbering();
- // Distribute the mesh
- ALE::Obj<partition_type> partition =
- DistributionType::distributeMeshV(origSieveMesh, newSieveMesh,
- renumbering,
- sendMeshOverlap, recvMeshOverlap);
- if (origSieveMesh->debug()) {
- std::cout << "["<<commRank<<"]: Mesh Renumbering:"
- << std::endl;
- for (SieveMesh::renumbering_type::const_iterator r_iter = renumbering.begin();
- r_iter != renumbering.end();
- ++r_iter) {
- std::cout << "["<<commRank<<"]: global point "
- << r_iter->first << " --> " << " local point "
- << r_iter->second << std::endl;
- } // for
- } // if
-
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Checking the overlap." << journal::endl;
- } // if
-
- // Check overlap
- int localSendOverlapSize = 0, sendOverlapSize;
- int localRecvOverlapSize = 0, recvOverlapSize;
- const int commSize = sendMeshOverlap->commSize();
- for (int p = 0; p < commSize; ++p) {
- localSendOverlapSize += sendMeshOverlap->getConeSize(p);
- localRecvOverlapSize += recvMeshOverlap->getSupportSize(p);
- } // for
- MPI_Allreduce(&localSendOverlapSize, &sendOverlapSize, 1, MPI_INT, MPI_SUM,
- sendMeshOverlap->comm());
- MPI_Allreduce(&localRecvOverlapSize, &recvOverlapSize, 1, MPI_INT, MPI_SUM,
- recvMeshOverlap->comm());
- if (sendOverlapSize != recvOverlapSize) {
- std::cout <<"["<<sendMeshOverlap->commRank()<<"]: Size mismatch " <<
- sendOverlapSize << " != " << recvOverlapSize << std::endl;
- sendMeshOverlap->view("Send Overlap");
- recvMeshOverlap->view("Recv Overlap");
- throw ALE::Exception("Invalid Overlap");
- } // if
-
- // Distribute the coordinates
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Distribution the vertex coordinates." << journal::endl;
- } // if
-
- const ALE::Obj<RealSection>& coordinates =
- origSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- const ALE::Obj<RealSection>& parallelCoordinates =
- newSieveMesh->getRealSection("coordinates");
- assert(!parallelCoordinates.isNull());
-
- newSieveMesh->setupCoordinates(parallelCoordinates);
- DistributionType::distributeSection(coordinates, partition, renumbering,
- sendMeshOverlap, recvMeshOverlap,
- parallelCoordinates);
-
- // Distribute other sections
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Distribution other sections." << journal::endl;
- } // if
-
- if (origSieveMesh->getRealSections()->size() > 1) {
- ALE::Obj<std::set<std::string> > names = origSieveMesh->getRealSections();
- assert(!names.isNull());
- int n = 0;
-
- const std::set<std::string>::const_iterator namesBegin = names->begin();
- const std::set<std::string>::const_iterator namesEnd = names->end();
- for (std::set<std::string>::const_iterator n_iter = namesBegin;
- n_iter != namesEnd;
- ++n_iter) {
- if (*n_iter == "coordinates") continue; // already copied
- if (*n_iter == "replaced_cells") continue; // ignore
- std::cerr << "ERROR: Did not distribute real section '" << *n_iter
- << "'." << std::endl;
- ++n;
- } // if
- if (n)
- throw std::logic_error("Need to distribute more real sections");
- }
-
- if (origSieveMesh->getIntSections()->size() > 0) {
- ALE::Obj<std::set<std::string> > names = origSieveMesh->getIntSections();
- const int numVertices = newSieveMesh->depthStratum(0)->size();
- const int numCells = newSieveMesh->heightStratum(0)->size();
- assert(!names.isNull());
-
- std::set<std::string>::const_iterator namesBegin = names->begin();
- std::set<std::string>::const_iterator namesEnd = names->end();
- for (std::set<std::string>::const_iterator n_iter = namesBegin;
- n_iter != namesEnd;
- ++n_iter) {
- const ALE::Obj<IntSection>& origSection =
- origSieveMesh->getIntSection(*n_iter);
- assert(!origSection.isNull());
- const ALE::Obj<IntSection>& newSection =
- newSieveMesh->getIntSection(*n_iter);
- assert(!newSection.isNull());
-
- // We assume all integer sections are supported on either cells or vertices
- SieveMesh::point_type firstPoint;
- IntSection::chart_type::const_iterator chartEnd = origSection->getChart().end();
- for(IntSection::chart_type::const_iterator c_iter = origSection->getChart().begin(); c_iter != chartEnd; ++c_iter) {
- if (origSection->getFiberDimension(*c_iter)) {
- firstPoint = *c_iter;
- break;
- } // if
- }
- if (origSieveMesh->height(firstPoint) == 0) {
- newSection->setChart(IntSection::chart_type(0, parallelCoordinates->getChart().min()));
- } else {
- newSection->setChart(parallelCoordinates->getChart());
- } // if/else
- DistributionType::distributeSection(origSection, partition, renumbering,
- sendMeshOverlap, recvMeshOverlap,
- newSection);
-#if 0 // DEBUGGING
- std::string serialName("Serial ");
- std::string parallelName("Parallel ");
- serialName += *n_iter;
- parallelName += *n_iter;
- origSection->view(serialName.c_str());
- newSection->view(parallelName.c_str());
-#endif
- } // for
- } // if
- if (origSieveMesh->getArrowSections()->size() > 1)
- throw std::logic_error("Need to distribute more arrow sections");
-
- // Distribute labels
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Distributing labels." << journal::endl;
- } // if
-
- const SieveMesh::labels_type& labels = origSieveMesh->getLabels();
- const SieveMesh::labels_type::const_iterator labelsBegin = labels.begin();
- const SieveMesh::labels_type::const_iterator labelsEnd = labels.end();
-
- for (SieveMesh::labels_type::const_iterator l_iter = labelsBegin;
- l_iter != labelsEnd;
- ++l_iter) {
- if (newSieveMesh->hasLabel(l_iter->first)) continue;
- const ALE::Obj<SieveMesh::label_type>& origLabel = l_iter->second;
- assert(!origLabel.isNull());
- const ALE::Obj<SieveMesh::label_type>& newLabel =
- newSieveMesh->createLabel(l_iter->first);
- assert(!newLabel.isNull());
-
-#ifdef IMESH_NEW_LABELS
- newLabel->setChart(newSieve->getChart());
- // Size the local mesh
- partitioner_type::sizeLocalSieveV(origLabel, partition, renumbering,
- newLabel);
- // Create the remote meshes
- DistributionType::completeConesV(origLabel, newLabel, renumbering,
- sendMeshOverlap, recvMeshOverlap);
- // Create the local mesh
- partitioner_type::createLocalSieveV(origLabel, partition, renumbering,
- newLabel);
- newLabel->symmetrize();
-#else
- DistributionType::distributeLabelV(newSieveMesh->getSieve(), origLabel, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newLabel);
-#if 0 // DEBUGGING
- std::string serialName("Serial ");
- std::string parallelName("Parallel ");
- serialName += l_iter->first;
- parallelName += l_iter->first;
- origLabel->view(serialName.c_str());
- newLabel->view(parallelName.c_str());
-#endif
-#endif
- } // for
-
- // Create the parallel overlap
- if (0 == commRank) {
- info << journal::at(__HERE__)
- << "Creating the parallel overlap." << journal::endl;
- } // if
-
- ALE::Obj<SieveMesh::send_overlap_type> sendParallelMeshOverlap =
- newSieveMesh->getSendOverlap();
- assert(!sendParallelMeshOverlap.isNull());
- ALE::Obj<SieveMesh::recv_overlap_type> recvParallelMeshOverlap =
- newSieveMesh->getRecvOverlap();
- assert(!recvParallelMeshOverlap.isNull());
-
- // Can I figure this out in a nicer way?
- ALE::SetFromMap<std::map<point_type,point_type> > globalPoints(renumbering);
- ALE::OverlapBuilder<>::constructOverlap(globalPoints, renumbering,
- sendParallelMeshOverlap,
- recvParallelMeshOverlap);
- newSieveMesh->setCalculatedOverlap(true);
-
-#if 0 // DEBUGGING
- sendParallelMeshOverlap->view("SEND OVERLAP");
- recvParallelMeshOverlap->view("RECV OVERLAP");
-#endif
-
- PYLITH_METHOD_END;
-} // distribute
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.hh 2013-05-11 03:43:34 UTC (rev 22042)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.hh 2013-05-11 04:08:39 UTC (rev 22043)
@@ -66,20 +66,6 @@
void write(meshio::DataWriter<topology::Mesh, topology::Field<topology::Mesh> >* const writer,
const topology::Mesh& mesh);
-// PRIVATE //////////////////////////////////////////////////////////////
-private :
-
- /** Distribute mesh among processors.
- *
- * @param newMesh Distributed mesh (result).
- * @param origMesh Mesh to distribute.
- */
- template<typename DistributionType>
- static
- void
- _distribute(topology::Mesh* const newMesh,
- const topology::Mesh& origMesh);
-
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc 2013-05-11 03:43:34 UTC (rev 22042)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc 2013-05-11 04:08:39 UTC (rev 22043)
@@ -211,51 +211,6 @@
const PylithScalar lengthScale = normalizer.lengthScale();
PetscErrorCode err;
-#if 1 // :TODO: REMOVE SIEVE STUFF
- // Get coordinates (currently dimensioned).
- if (!_mesh.isNull()) {
- const ALE::Obj<RealSection>& coordsSection = _mesh->getRealSection("coordinates");
- assert(!coordsSection.isNull());
-
- // Get field for dimensioned coordinates.
- const ALE::Obj<RealSection>& coordsDimSection =
- _mesh->getRealSection("coordinates_dimensioned");
- assert(!coordsDimSection.isNull());
- coordsDimSection->setAtlas(coordsSection->getAtlas());
- coordsDimSection->allocateStorage();
- coordsDimSection->setBC(coordsSection->getBC());
-
- const ALE::Obj<SieveMesh::label_sequence>& vertices = _mesh->depthStratum(0);
- assert(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
-
- PylithScalar coordsVertex[3];
- for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- const int spaceDim = coordsSection->getFiberDimension(*v_iter);
- assert(spaceDim <= 3);
- const PylithScalar* coordsDimVertex = coordsSection->restrictPoint(*v_iter);
-
- // Update section with dimensioned coordinates
- assert(spaceDim == coordsDimSection->getFiberDimension(*v_iter));
- coordsDimSection->updatePoint(*v_iter, coordsDimVertex);
-
- // Copy coordinates to array for nondimensionalization.
- for (int i=0; i < spaceDim; ++i)
- coordsVertex[i] = coordsDimVertex[i];
-
- // Nondimensionalize original coordinates.
- normalizer.nondimensionalize(&coordsVertex[0], spaceDim, lengthScale);
-
- // Update section with nondimensional coordinates
- assert(spaceDim == coordsSection->getFiberDimension(*v_iter));
- coordsSection->updatePoint(*v_iter, coordsVertex);
- } // for
- }
-#endif
-
assert(_newMesh);
err = DMGetCoordinatesLocal(_newMesh, &coordVec);PYLITH_CHECK_ERROR(err);assert(coordVec);
// There does not seem to be an advantage to calling nondimensionalize()
More information about the CIG-COMMITS
mailing list