[cig-commits] r20894 - in short/3D/PyLith/trunk/libsrc/pylith: meshio topology

[knepley at geodynamics.org]

Author: knepley
Date: 2012-10-23 21:08:16 -0700 (Tue, 23 Oct 2012)
New Revision: 20894

Modified:
   short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc
Log:
Now non-dimensionalize new coordinates

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc	2012-10-24 02:27:33 UTC (rev 20893)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc	2012-10-24 04:08:16 UTC (rev 20894)
@@ -119,10 +119,12 @@
     if (dmMesh) {
       PetscSection coordSection;
       Vec          coordinates;
+      PetscReal    lengthScale;
       PetscInt     vStart, vEnd, vMax, verticesSize, dim, dimLocal = 0;
 
       const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
       assert(cs);
+      err = DMComplexGetScale(dmMesh, PETSC_UNIT_LENGTH, &lengthScale);CHECK_PETSC_ERROR(err);
       err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
       err = DMGetCoordinatesLocal(dmMesh, &coordinates);CHECK_PETSC_ERROR(err);
       err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
@@ -151,6 +153,7 @@
       }
       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);
       err = PetscViewerHDF5PushGroup(_viewer, "/geometry");CHECK_PETSC_ERROR(err);
       err = VecView(coordVec, _viewer);CHECK_PETSC_ERROR(err);

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc	2012-10-24 02:27:33 UTC (rev 20893)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc	2012-10-24 04:08:16 UTC (rev 20894)
@@ -142,17 +142,19 @@
     PetscVec coordinatesVector = coordinates.vector(context);
     assert(coordinatesVector);
 #else
-    DM       coordDM;
-    Vec      coordVec, coordinatesVector;
-    PetscInt globalSize;
+    DM        coordDM;
+    Vec       coordVec, coordinatesVector;
+    PetscReal lengthScale;
+    PetscInt  globalSize;
 
-    /* Should use the coordinate DM here (once I put it in) */
+    err = DMComplexGetScale(dmMesh, PETSC_UNIT_LENGTH, &lengthScale);CHECK_PETSC_ERROR(err);
     err = DMGetCoordinateDM(dmMesh, &coordDM);CHECK_PETSC_ERROR(err);
     err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
     err = DMGetGlobalVector(coordDM, &coordinatesVector);CHECK_PETSC_ERROR(err);
     err = VecGetSize(coordinatesVector, &globalSize);CHECK_PETSC_ERROR(err);
     err = DMLocalToGlobalBegin(coordDM, coordVec, INSERT_VALUES, coordinatesVector);CHECK_PETSC_ERROR(err);
     err = DMLocalToGlobalEnd(coordDM, coordVec, INSERT_VALUES, coordinatesVector);CHECK_PETSC_ERROR(err);
+    err = VecScale(coordinatesVector, lengthScale);CHECK_PETSC_ERROR(err);
 #endif
 
     const std::string& filenameVertices = _datasetFilename("vertices");

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc	2012-10-24 02:27:33 UTC (rev 20893)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.cc	2012-10-24 04:08:16 UTC (rev 20894)
@@ -116,14 +116,19 @@
 void 
 pylith::topology::Mesh::nondimensionalize(const spatialdata::units::Nondimensional& normalizer)
 { // initialize
-  //TODO Change DM coordinates here, store lengthScale for output routines
 
   // Get coordinates (currently dimensioned).
   assert(!_mesh.isNull());
   const ALE::Obj<RealSection>& coordsSection =
     _mesh->getRealSection("coordinates");
   assert(!coordsSection.isNull());
+  Vec coordVec, coordDimVec;
+  PetscErrorCode err;
 
+  assert(_newMesh);
+  err = DMGetCoordinatesLocal(_newMesh, &coordVec);CHECK_PETSC_ERROR(err);
+  assert(coordVec);
+
   // Get field for dimensioned coordinates.
   const ALE::Obj<RealSection>& coordsDimSection =
     _mesh->getRealSection("coordinates_dimensioned");
@@ -141,6 +146,10 @@
   const SieveMesh::label_sequence::iterator verticesEnd = 
     vertices->end();
 
+  // There does not seem to be an advantage to calling nondimensionalize()
+  err = VecScale(coordVec, 1.0/lengthScale);CHECK_PETSC_ERROR(err);
+  err = DMComplexSetScale(_newMesh, PETSC_UNIT_LENGTH, lengthScale);CHECK_PETSC_ERROR(err);
+
   PylithScalar coordsVertex[3];
   for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
       v_iter != verticesEnd;