[cig-commits] r8395 - in short/3D/PyLith/trunk: . libsrc/faults

[knepley at geodynamics.org]

Author: knepley
Date: 2007-12-05 14:53:10 -0800 (Wed, 05 Dec 2007)
New Revision: 8395

Modified:
   short/3D/PyLith/trunk/TODO
   short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
Log:
More debugging


Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO	2007-12-05 22:00:33 UTC (rev 8394)
+++ short/3D/PyLith/trunk/TODO	2007-12-05 22:53:10 UTC (rev 8395)
@@ -278,6 +278,8 @@
 
 5. Throw out cohesive midnodes from censored depth (not easy since they are mixed in with new vertices, maybe change numbering)
 
+   I think you can do this by placing them after cohesive cells
+
 6. Coarsening for PyLith meshes
 
 7. Time fixes
@@ -294,7 +296,7 @@
 
 10. Move events to Python (make wrappers for events)
 
-11. For Schur complement directly for mass matrix and Lagrang multiplier problem
+11. For Schur complement directly for mass matrix and Lagrange multiplier problem
 
 ======================================================================
 QUESTIONS FOR LEIF

Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc	2007-12-05 22:00:33 UTC (rev 8394)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc	2007-12-05 22:53:10 UTC (rev 8395)
@@ -515,26 +515,47 @@
   }
 
   // More checking
+  const bool                         firstFault    = !mesh->hasRealSection("replacedCells");
   const ALE::Obj<real_section_type>& replacedCells = mesh->getRealSection("replacedCells");
   PointSet cellNeighbors;
 
-  replacedCells->setFiberDimension(mesh->heightStratum(0), 1);
-  replacedCells->allocatePoint();
+  if (firstFault) {
+    replacedCells->setFiberDimension(mesh->heightStratum(0), 1);
+    replacedCells->allocatePoint();
+  }
   for(PointSet::const_iterator c_iter = noReplaceCells.begin(); c_iter != noReplaceCells.end(); ++c_iter) {
     const double minusOne = -1.0;
 
-    replacedCells->updatePoint(*c_iter, &minusOne);
+    if (replacedCells->restrictPoint(*c_iter)[0] == 0.0) {
+      replacedCells->updatePoint(*c_iter, &minusOne);
+    } else {
+      const double minusTwo = -2.0;
+
+      replacedCells->updatePoint(*c_iter, &minusTwo);
+    }
   }
   for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
     if (replaceCellsBase.find(*c_iter) != replaceCellsBase.end()) {
       const double one = 1.0;
 
-      replacedCells->updatePoint(*c_iter, &one);
+      if (replacedCells->restrictPoint(*c_iter)[0] == 0.0) {
+        replacedCells->updatePoint(*c_iter, &one);
+      } else {
+        const double two = 2.0;
+
+        replacedCells->updatePoint(*c_iter, &two);
+      }
       continue;
     }
     const double ten = 10.0;
 
-    replacedCells->updatePoint(*c_iter, &ten);
+    if (replacedCells->restrictPoint(*c_iter)[0] == 0.0) {
+      replacedCells->updatePoint(*c_iter, &ten);
+    } else {
+        const double twenty = 20.0;
+
+        replacedCells->updatePoint(*c_iter, &twenty);
+    }
     // There should be a way to check for boundary elements
     if (mesh->getDimension() == 1) {
       if (cellNeighbors.size() > 2) {