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

[knepley at geodynamics.org]

Author: knepley
Date: 2012-01-26 20:26:11 -0800 (Thu, 26 Jan 2012)
New Revision: 19487

Modified:
   short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc
Log:
Test code for Jacobian (disabled)


Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc	2012-01-27 01:13:29 UTC (rev 19486)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc	2012-01-27 04:26:11 UTC (rev 19487)
@@ -126,7 +126,7 @@
     CHECK_PETSC_ERROR(err);
 
     const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
-    
+
     err = VTKViewer::writeHeader(sieveMesh, _viewer);
     CHECK_PETSC_ERROR(err);
     //std::cout << "Wrote header for " << filename << std::endl;

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc	2012-01-27 01:13:29 UTC (rev 19486)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc	2012-01-27 04:26:11 UTC (rev 19487)
@@ -128,6 +128,32 @@
     CHECK_PETSC_ERROR(err);
     err = MatAssemblyEnd(_matrix, MAT_FINAL_ASSEMBLY);
     CHECK_PETSC_ERROR(err);
+#if 0
+    // Check for empty row
+    const PetscInt *cols;
+    PetscInt        rStart, rEnd, ncols;
+
+    err = MatGetOwnershipRange(_matrix, &rStart, &rEnd);CHECK_PETSC_ERROR(err);
+    for(PetscInt r = rStart; r < rEnd; ++r) {
+      PetscInt c;
+
+      err = MatGetRow(_matrix,r, &ncols, &cols, PETSC_NULL);CHECK_PETSC_ERROR(err);
+      if (!ncols) {
+        std::ostringstream msg;
+        msg << "ERROR: Empty row " << r << " in ["<<rStart<<","<<rEnd<<")" << std::endl;
+        throw std::runtime_error(msg.str().c_str());
+      }
+      for(c = 0; c < ncols; ++c) {
+        if (cols[c] == r) break;
+      }
+      if (c == ncols) {
+        std::ostringstream msg;
+        msg << "ERROR: Row " << r << " in ["<<rStart<<","<<rEnd<<") is missing diagonal element" << std::endl;
+        throw std::runtime_error(msg.str().c_str());
+      }
+      err = MatRestoreRow(_matrix,r, &ncols, &cols, PETSC_NULL);CHECK_PETSC_ERROR(err);
+    }
+#endif
   } else if (0 == strcmp(mode, "flush_assembly")) {
     err = MatAssemblyBegin(_matrix, MAT_FLUSH_ASSEMBLY);
     CHECK_PETSC_ERROR(err);