[cig-commits] r19891 - in short/3D/PyLith/branches/pylith-scecdynrup: examples/3d/hex8 libsrc/pylith/meshio libsrc/pylith/topology
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 28 16:09:46 PDT 2012
Author: brad
Date: 2012-03-28 16:09:46 -0700 (Wed, 28 Mar 2012)
New Revision: 19891
Modified:
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/MeshBuilder.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/MeshBuilder.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputManager.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Field.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Mesh.hh
Log:
Merge from trunk.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt 2012-03-28 23:06:42 UTC (rev 19890)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt 2012-03-28 23:09:46 UTC (rev 19891)
@@ -7,7 +7,7 @@
# coordsys component of the OutputSolnPoints object. The points will
# be transformed into the coordinate system of the mesh before
# interpolation.
-0.0001 0.0001 -500.0
-0.0001 0.0001 -1500.0
-0.0001 0.0001 -2500.0
-0.0001 0.0001 -3500.0
+0.0 0.0 -500.0
+0.0 0.0 -1500.0
+0.0 0.0 -2500.0
+0.0 0.0 -3500.0
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/MeshBuilder.cc 2012-03-28 23:06:42 UTC (rev 19890)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/MeshBuilder.cc 2012-03-28 23:09:46 UTC (rev 19891)
@@ -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/branches/pylith-scecdynrup/libsrc/pylith/meshio/MeshBuilder.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/MeshBuilder.hh 2012-03-28 23:06:42 UTC (rev 19890)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/MeshBuilder.hh 2012-03-28 23:09:46 UTC (rev 19891)
@@ -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/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputManager.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputManager.cc 2012-03-28 23:06:42 UTC (rev 19890)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputManager.cc 2012-03-28 23:09:46 UTC (rev 19891)
@@ -248,7 +248,7 @@
throw std::logic_error("Unknown field type");
} // switch
- if (0 == _fields)
+ if (!_fields)
_fields = new topology::Fields<field_type>(fieldIn.mesh());
if (!_fields->hasField(fieldName.c_str())) {
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-28 23:06:42 UTC (rev 19890)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-28 23:09:46 UTC (rev 19891)
@@ -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,22 +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;
- 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;
@@ -138,27 +124,62 @@
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);
+ err = DMMeshInterpolationAddPoints(dm, numPoints, (PetscReal*)&pointsNondim[0], _interpolator);CHECK_PETSC_ERROR(err);
+ const PetscBool pointsAllProcs = PETSC_TRUE;
+ err = DMMeshInterpolationSetUp(dm, _interpolator, pointsAllProcs);CHECK_PETSC_ERROR(err);
+ err = DMDestroy(&dm);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());
- assert(_pointsMesh->sieveMesh()->hasRealSection("coordinates"));
- const ALE::Obj<topology::Mesh::RealSection>& coordinatesSection = _pointsMesh->sieveMesh()->getRealSection("coordinates");
- assert(!coordinatesSection.isNull());
- const PylithScalar* coordinates = coordinatesSection->restrictSpace();
- err = DMMeshInterpolationAddPoints(dm, numPoints, (PetscReal*)coordinates,
- _interpolator);CHECK_PETSC_ERROR(err);
- err = DMMeshInterpolationSetUp(dm, _interpolator);CHECK_PETSC_ERROR(err);
- err = DMDestroy(&dm);CHECK_PETSC_ERROR(err);
- CHECK_PETSC_ERROR(err);
+ 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
@@ -234,7 +255,7 @@
topology::Field<topology::Mesh>& fieldInterp =
_fields->get("buffer (interpolated)");
- if (vertices->size()*fiberDim != fieldInterp.sectionSize()) {
+ if (fieldInterp.section().isNull() || vertices->size()*fiberDim != fieldInterp.sectionSize()) {
fieldInterp.newSection(vertices, fiberDim);
fieldInterp.label(field.label());
fieldInterp.allocate();
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Field.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Field.hh 2012-03-28 23:06:42 UTC (rev 19890)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Field.hh 2012-03-28 23:09:46 UTC (rev 19891)
@@ -414,6 +414,8 @@
Metadata _metadata;
const mesh_type& _mesh; ///< Mesh associated with section.
ALE::Obj<section_type> _section; ///< Real section with data.
+ PetscSection _newSection;
+ Vec _newLocalVec;
scatter_map_type _scatters; ///< Collection of scatters.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Mesh.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Mesh.hh 2012-03-28 23:06:42 UTC (rev 19890)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Mesh.hh 2012-03-28 23:09:46 UTC (rev 19891)
@@ -197,6 +197,7 @@
private :
ALE::Obj<SieveMesh> _mesh; ///< Sieve mesh.
+ DM _newMesh;
spatialdata::geocoords::CoordSys* _coordsys; ///< Coordinate system.
MPI_Comm _comm; ///< MPI communicator for mesh.
bool _debug; ///< Debugging flag for mesh.
More information about the CIG-COMMITS
mailing list