[cig-commits] r19842 - short/3D/PyLith/trunk/libsrc/pylith/meshio
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 21 15:49:46 PDT 2012
Author: brad
Date: 2012-03-21 15:49:45 -0700 (Wed, 21 Mar 2012)
New Revision: 19842
Modified:
short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.hh
short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc
Log:
Reworked setup of point interpolation so point mesh is parallel.
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc 2012-03-21 22:23:42 UTC (rev 19841)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc 2012-03-21 22:49:45 UTC (rev 19842)
@@ -45,11 +45,12 @@
const int numCells,
const int numCorners,
const int meshDim,
- const bool interpolate)
+ const bool interpolate,
+ const bool isParallel)
{ // buildMesh
- assert(0 != mesh);
+ assert(mesh);
- assert(0 != coordinates);
+ assert(coordinates);
MPI_Comm comm = mesh->comm();
int dim = meshDim;
int rank = 0;
@@ -87,7 +88,7 @@
//logger.setDebug(1);
logger.stagePush("MeshCreation");
- if (0 == rank) {
+ if (0 == rank || isParallel) {
assert(coordinates->size() == numVertices*spaceDim);
assert(cells.size() == numCells*numCorners);
if (!interpolate) {
@@ -149,11 +150,11 @@
for(int c = 0; c < numCells+numVertices; ++c) {
height->setConeSize(c, 1);
depth->setConeSize(c, 1);
- }
- if (numCells+numVertices)
+ } // for
+ if (numCells+numVertices) {
height->setSupportSize(0, numCells+numVertices);
- if (numCells+numVertices)
depth->setSupportSize(0, numCells+numVertices);
+ } // if
height->allocate();
depth->allocate();
#endif
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.hh 2012-03-21 22:23:42 UTC (rev 19841)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.hh 2012-03-21 22:49:45 UTC (rev 19842)
@@ -60,6 +60,9 @@
* @param numCells Number of cells.
* @param numCorners Number of vertices per cell.
* @param meshDim Dimension of cells in mesh.
+ * @param interpolate Create interpolated mesh.
+ * @param isParallel Create parallel mesh if true, otherwise only build
+ * mesh on proc 0.
*/
static
void buildMesh(topology::Mesh* mesh,
@@ -70,7 +73,8 @@
const int numCells,
const int numCorners,
const int meshDim,
- const bool interpolate);
+ const bool interpolate,
+ const bool isParallel =false);
/** Build fault mesh topology.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-21 22:23:42 UTC (rev 19841)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-21 22:49:45 UTC (rev 19842)
@@ -90,20 +90,19 @@
assert(mesh);
assert(points);
+ assert(!_interpolator); // Insure clean starting point
+
_mesh = mesh;
- // Create mesh without cells for points.
- const int meshDim = 0;
- delete _pointsMesh; _pointsMesh = new topology::Mesh(meshDim);
- _pointsMesh->createSieveMesh(meshDim);
- assert(_pointsMesh);
+ // Create nondimensionalized array of point locations
+ const int size = numPoints * spaceDim;
+ scalar_array pointsNondim(size);
+ for (int i=0; i < size; ++i) {
+ pointsNondim[i] = points[i] / normalizer.lengthScale();
+ } // for
- const spatialdata::geocoords::CoordSys* csMesh = mesh->coordsys();
- assert(csMesh);
- assert(csMesh->spaceDim() == spaceDim);
-
#if 1 // DEBUGGING
- std::cout << "OUTPUT SOLN POINTS" << std::endl;
+ std::cout << "OUTPUT SOLN POINTS (dimensioned)" << std::endl;
for (int i=0; i < numPoints; ++i) {
for (int iDim=0; iDim < spaceDim; ++iDim) {
std::cout << " " << points[i*spaceDim+iDim];
@@ -112,23 +111,9 @@
} // for
#endif
- scalar_array pointsArray(points, numPoints*spaceDim);
- int_array cells(numPoints);
- for (int i=0; i < numPoints; ++i)
- cells[i] = i;
- const int numCells = numPoints;
- const int numCorners = 1;
- const bool interpolate = false;
- // Build mesh creates mesh on proc 0
- MeshBuilder::buildMesh(_pointsMesh,
- &pointsArray, numPoints, spaceDim,
- cells, numCells, numCorners, meshDim,
- interpolate);
- _pointsMesh->coordsys(_mesh->coordsys());
- _pointsMesh->nondimensionalize(normalizer);
-#if 1 // DEBUGGING
- _pointsMesh->view("POINTS MESH");
-#endif
+ const spatialdata::geocoords::CoordSys* csMesh = mesh->coordsys();
+ assert(csMesh);
+ assert(csMesh->spaceDim() == spaceDim);
// Setup interpolator object
DM dm;
@@ -139,29 +124,66 @@
err = DMMeshSetMesh(dm, _mesh->sieveMesh());CHECK_PETSC_ERROR(err);
err = DMMeshInterpolationCreate(dm, &_interpolator);CHECK_PETSC_ERROR(err);
- const spatialdata::geocoords::CoordSys* cs = _pointsMesh->coordsys();
- assert(cs);
- assert(cs->spaceDim() == spaceDim);
+ err = DMMeshInterpolationSetDim(dm, spaceDim, _interpolator);CHECK_PETSC_ERROR(err);
- err = DMMeshInterpolationSetDim(dm, spaceDim,
- _interpolator);CHECK_PETSC_ERROR(err);
-
- assert(!_pointsMesh->sieveMesh().isNull());
- assert(_pointsMesh->sieveMesh()->hasRealSection("coordinates"));
- const ALE::Obj<topology::Mesh::RealSection>& coordinatesSection = _pointsMesh->sieveMesh()->getRealSection("coordinates");
- assert(!coordinatesSection.isNull());
- assert(0 == coordinatesSection->sizeWithBC() % spaceDim);
- const int numPointsLocal = coordinatesSection->sizeWithBC() / spaceDim;
- const PylithScalar* coordinates = (numPointsLocal) ? coordinatesSection->restrictSpace() : 0;
- err = DMMeshInterpolationAddPoints(dm, numPointsLocal, (PetscReal*)coordinates,
- _interpolator);CHECK_PETSC_ERROR(err);
+ err = DMMeshInterpolationAddPoints(dm, numPoints, (PetscReal*)&pointsNondim[0], _interpolator);CHECK_PETSC_ERROR(err);
+ const bool pointsAllProcs = true;
+#if 0 // WAITING FOR MATT TO IMPLEMENT
+ err = DMMeshInterpolationSetUp(dm, _interpolator, pointsAllProcs);CHECK_PETSC_ERROR(err);
+#else
err = DMMeshInterpolationSetUp(dm, _interpolator);CHECK_PETSC_ERROR(err);
+#endif
err = DMDestroy(&dm);CHECK_PETSC_ERROR(err);
- CHECK_PETSC_ERROR(err);
+ // Create mesh corresponding to points.
+ const int meshDim = 0;
+ delete _pointsMesh; _pointsMesh = new topology::Mesh(meshDim);
+ _pointsMesh->createSieveMesh(meshDim);
+ assert(_pointsMesh);
+
+ const int numPointsLocal = _interpolator->n;
+ PylithScalar* pointsLocal = 0;
+ err = VecGetArray(_interpolator->coords, &pointsLocal);CHECK_PETSC_ERROR(err);
+ scalar_array pointsArray(numPointsLocal*spaceDim);
+ const int sizeLocal = numPointsLocal*spaceDim;
+ for (int i=0; i < sizeLocal; ++i) {
+ // Must scale by length scale because we gave interpolator nondimensioned coordinates
+ pointsArray[i] = pointsLocal[i]*normalizer.lengthScale();
+ } // for
+ int_array cells(numPointsLocal);
+ for (int i=0; i < numPointsLocal; ++i)
+ cells[i] = i;
+ const int numCells = numPointsLocal;
+ const int numCorners = 1;
+ const bool interpolate = false;
+ const bool isParallel = true;
+ MeshBuilder::buildMesh(_pointsMesh, &pointsArray, numPointsLocal, spaceDim,
+ cells, numCells, numCorners, meshDim, interpolate, isParallel);
+ err = VecRestoreArray(_interpolator->coords, &pointsLocal);CHECK_PETSC_ERROR(err);
+
+ // Set coordinate system and create nondimensionalized coordinates
+ _pointsMesh->coordsys(_mesh->coordsys());
+ _pointsMesh->nondimensionalize(normalizer);
+
+ // Create empty overlaps. MATT - IS THIS OKAY?
+ assert(!_pointsMesh->sieveMesh().isNull());
+ ALE::Obj<SieveSubMesh::send_overlap_type> sendOverlap = _pointsMesh->sieveMesh()->getSendOverlap();
+ assert(!sendOverlap.isNull());
+ ALE::Obj<SieveSubMesh::recv_overlap_type> recvOverlap = _pointsMesh->sieveMesh()->getRecvOverlap();
+ assert(!recvOverlap.isNull());
+ SieveMesh::renumbering_type emptyRenumbering;
+ ALE::SetFromMap<SieveMesh::renumbering_type> emptyPoints(emptyRenumbering);
+ ALE::OverlapBuilder<>::constructOverlap(emptyPoints, emptyRenumbering, sendOverlap, recvOverlap);
+ _pointsMesh->sieveMesh()->setCalculatedOverlap(true);
+
+#if 1 // DEBUGGING
+ _pointsMesh->view("POINTS MESH");
+#endif
+
if (!_fields) {
_fields = new topology::Fields<topology::Field<topology::Mesh> >(*_pointsMesh);
} // if
+
} // setupInterpolator
More information about the CIG-COMMITS
mailing list