[cig-commits] r19142 - in short/3D/PyLith/trunk: examples/3d/hex8 libsrc/pylith/meshio libsrc/pylith/topology

[knepley at geodynamics.org]

Author: knepley
Date: 2011-11-02 06:59:56 -0700 (Wed, 02 Nov 2011)
New Revision: 19142

Modified:
   short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
   short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/MeshRefiner.cc
Log:
Added reordering to hex8, Fixed DM wrapper for interpolation, added debugging for refinement


Modified: short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg	2011-11-02 01:48:33 UTC (rev 19141)
+++ short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg	2011-11-02 13:59:56 UTC (rev 19142)
@@ -30,6 +30,10 @@
 # Change the default mesh reader to the CUBIT reader.
 reader = pylith.meshio.MeshIOCubit
 
+# Optimize ordering of mesh cells and vertices using reverse
+# Cuthill-KcKee algorithm.
+reorder_mesh = True
+
 [pylithapp.mesh_generator.reader]
 # Set filename of mesh to import.
 filename = mesh/box_hex8_1000m.exo

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc	2011-11-02 01:48:33 UTC (rev 19141)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc	2011-11-02 13:59:56 UTC (rev 19142)
@@ -52,15 +52,17 @@
 { // deallocate
   OutputManager<topology::Mesh, topology::Field<topology::Mesh> >::deallocate();
 
-#if 0 // :MATT: Need DM wrapper for mesh
   if (_interpolator) {
     assert(_mesh);
     const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
-    PetscErrorCode err = 
-      DMMeshInterpolationDestroy(sieveMesh, &_interpolator);
-    CHECK_PETSC_ERROR(err);
+    DM dm;
+    PetscErrorCode err = 0;
+
+    err = DMMeshCreate(sieveMesh->comm(), &dm);CHECK_PETSC_ERROR(err);
+    err = DMMeshSetMesh(dm, sieveMesh);CHECK_PETSC_ERROR(err);
+    err = DMMeshInterpolationDestroy(dm, &_interpolator);CHECK_PETSC_ERROR(err);
+    err = DMDestroy(&dm);CHECK_PETSC_ERROR(err);
   } // if
-#endif
 
   _mesh = 0; // :TODO: Use shared pointer
   delete _pointsMesh; _pointsMesh = 0;
@@ -108,26 +110,23 @@
 			 interpolate);
   _pointsMesh->coordsys(_mesh->coordsys());
 
-#if 0 // :MATT: Need DM wrapper for mesh
   // Setup interpolator object
+  DM dm;
   PetscErrorCode err = 0;
 
-  err = DMMeshInterpolationCreate(_mesh->sieveMesh(), &_interpolator); 
-  CHECK_PETSC_ERROR(err);
+  err = DMMeshCreate(_mesh->sieveMesh()->comm(), &dm);CHECK_PETSC_ERROR(err);
+  err = DMMeshSetMesh(dm, _mesh->sieveMesh());CHECK_PETSC_ERROR(err);
+  err = DMMeshInterpolationCreate(dm, &_interpolator);CHECK_PETSC_ERROR(err);
   
-  const spatialdata::geocoords::CoordSys* cs = _pointsMesh().coordsys();
+  const spatialdata::geocoords::CoordSys* cs = _pointsMesh->coordsys();
   assert(0 != cs);
   assert(cs->spaceDim() == spaceDim);
 
-  err = DMMeshInterpolationSetDim(_mesh->sieveMesh(), spaceDim,
-				  _interpolator); CHECK_PETSC_ERROR(err);
-
-  err = DMMeshInterpolationAddPoints(_mesh->sieveMesh(), numPoints, points,
-				     _interpolator); CHECK_PETSC_ERROR(err);
-  
-  err = DMMeshInterpolationSetUp(_mesh->sieveMesh(), _interpolator);
+  err = DMMeshInterpolationSetDim(dm, spaceDim, _interpolator);CHECK_PETSC_ERROR(err);
+  err = DMMeshInterpolationAddPoints(dm, numPoints, (PetscReal *) points, _interpolator);CHECK_PETSC_ERROR(err);
+  err = DMMeshInterpolationSetUp(dm, _interpolator);CHECK_PETSC_ERROR(err);
+  err = DMDestroy(&dm);CHECK_PETSC_ERROR(err);
   CHECK_PETSC_ERROR(err);
-#endif
 } // createPointsMesh
 
 
@@ -176,16 +175,22 @@
   const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
   assert(!sieveMesh.isNull());
 
-#if 0 // :MATT: Need DM wrapper for mesh
+  DM dm;
+  SectionReal section;
   PetscErrorCode err = 0;
   
-  err = DMMeshInterpolationSetDof(sieveMesh, fiberDim, _interpolator);
-  CHECK_PETSC_ERROR(err);
-  
-  err = DMMeshInterpolationEvaluate(sieveMesh, field.section(), fieldInterpVec,
-				    _interpolator); CHECK_PETSC_ERROR(err);
-#endif  
+  err = DMMeshCreate(sieveMesh->comm(), &dm);CHECK_PETSC_ERROR(err);
+  err = DMMeshSetMesh(dm, sieveMesh);CHECK_PETSC_ERROR(err);
+  err = DMMeshInterpolationSetDof(dm, fiberDim, _interpolator);CHECK_PETSC_ERROR(err);
 
+  err = SectionRealCreate(sieveMesh->comm(), &section);CHECK_PETSC_ERROR(err);
+  err = SectionRealSetSection(section, field.section());CHECK_PETSC_ERROR(err);
+  err = SectionRealSetBundle(section, sieveMesh);CHECK_PETSC_ERROR(err);
+
+  err = DMMeshInterpolationEvaluate(dm, section, fieldInterpVec, _interpolator);CHECK_PETSC_ERROR(err);
+  err = SectionRealDestroy(&section);CHECK_PETSC_ERROR(err);
+  err = DMDestroy(&dm);CHECK_PETSC_ERROR(err);
+
   OutputManager<topology::Mesh, topology::Field<topology::Mesh> >::appendVertexField(t, fieldInterp, mesh);
 } // appendVertexField
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/MeshRefiner.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/MeshRefiner.cc	2011-11-02 01:48:33 UTC (rev 19141)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/MeshRefiner.cc	2011-11-02 13:59:56 UTC (rev 19142)
@@ -799,6 +799,12 @@
   // We have to do flexible assembly since we add the new vertices separately
   newSendOverlap->assemble();
   newRecvOverlap->assemble();
+  if (mesh->debug()) {
+    sendOverlap->view("OLD SEND OVERLAP");
+    recvOverlap->view("OLD RECV OVERLAP");
+    newSendOverlap->view("NEW SEND OVERLAP");
+    newRecvOverlap->view("NEW RECV OVERLAP");
+  }
 } // _calcNewOverlap