[cig-commits] r8035 - in short/3D/PyLith/trunk: libsrc/faults unittests/libtests/faults unittests/libtests/faults/data

knepley at geodynamics.org knepley at geodynamics.org
Wed Sep 26 20:32:37 PDT 2007


Author: knepley
Date: 2007-09-26 20:32:37 -0700 (Wed, 26 Sep 2007)
New Revision: 8035

Modified:
   short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataQuad4h.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataTri3f.cc
Log:
New algorithm for cohesive cells


Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc	2007-09-26 20:12:18 UTC (rev 8034)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc	2007-09-27 03:32:37 UTC (rev 8035)
@@ -39,6 +39,7 @@
   *fault = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
   const ALE::Obj<sieve_type> faultSieve = new sieve_type(sieve->comm(), 
                                                          sieve->debug());
+  const int  depth      = mesh->depth();
   const int  numCells   = mesh->heightStratum(0)->size();
   int        numCorners = 0;    // The number of vertices in a mesh cell
   int        faceSize   = 0;    // The number of vertices in a mesh face
@@ -49,7 +50,7 @@
   PointArray neighborVertices;
 
   if (!(*fault)->commRank()) {
-    numCorners = sieve->nCone(*mesh->heightStratum(0)->begin(), mesh->depth())->size();
+    numCorners = sieve->nCone(*mesh->heightStratum(0)->begin(), depth)->size();
     faceSize   = selection::numFaceVertices(*mesh->heightStratum(0)->begin(), mesh);
     indices    = new int[faceSize];
   }
@@ -73,7 +74,7 @@
   // Create a sieve which captures the fault
   for(PointSet::const_iterator fv_iter = fvBegin; fv_iter != fvEnd; ++fv_iter) {
     const ALE::Obj<sieveAlg::supportArray>& cells =
-      sieveAlg::nSupport(mesh, *fv_iter, mesh->depth());
+      sieveAlg::nSupport(mesh, *fv_iter, depth);
     const sieveAlg::supportArray::iterator cBegin = cells->begin();
     const sieveAlg::supportArray::iterator cEnd   = cells->end();
     
@@ -415,10 +416,91 @@
     }
     mesh->setValue(material, newPoint, materialId);
   } // for
-  // Replace all cells with a vertex on the fault that share a face with this one, or with one that does
+  // Replace all cells on a given side of the fault with a vertex on the fault
+#if 1
+  PointSet vReplaceCells;
+  PointSet vNoReplaceCells;
+
   for(PointSet::const_iterator v_iter = replaceVertices.begin(); v_iter != replaceVertices.end(); ++v_iter) {
+    const ALE::Obj<sieve_type::traits::supportSequence>& neighbors = sieve->support(*v_iter);
+    const sieve_type::traits::supportSequence::iterator  begin     = neighbors->begin();
+    const sieve_type::traits::supportSequence::iterator  end       = neighbors->end();
+    bool                                                 modified  = true;
+    int                                                  classifyTotal = neighbors->size();
+    int                                                  classifySize  = 0;
+
+    if (debug) {std::cout << "Checking fault vertex " << *v_iter << std::endl;}
+    vReplaceCells.clear();
+    vNoReplaceCells.clear();
+    for(sieve_type::traits::supportSequence::iterator n_iter = begin; n_iter != end; ++n_iter) {
+      if (replaceCells.find(*n_iter)   != replaceCells.end())   vReplaceCells.insert(*n_iter);
+      if (noReplaceCells.find(*n_iter) != noReplaceCells.end()) vNoReplaceCells.insert(*n_iter);
+      if (*n_iter >= firstCohesiveCell) classifyTotal--;
+    }
+    classifySize = vReplaceCells.size() + vNoReplaceCells.size();
+
+    while(modified && (classifySize < classifyTotal)) {
+      modified = false;
+      for(sieve_type::traits::supportSequence::iterator n_iter = begin; n_iter != end; ++n_iter) {
+        bool classified = false;
+
+        if (debug) {std::cout << "Checking neighbor " << *n_iter << std::endl;}
+        if (vReplaceCells.find(*n_iter)   != vReplaceCells.end()) {
+          if (debug) {std::cout << "  already in replaceCells" << std::endl;}
+          continue;
+        }
+        if (vNoReplaceCells.find(*n_iter) != vNoReplaceCells.end()) {
+          if (debug) {std::cout << "  already in noReplaceCells" << std::endl;}
+          continue;
+        }
+        if (*n_iter >= firstCohesiveCell) {
+          if (debug) {std::cout << "  already a cohesive cell" << std::endl;}
+          continue;
+        }
+        // If neighbor shares a face with anyone in replaceCells, then add
+        for(PointSet::const_iterator c_iter = vReplaceCells.begin(); c_iter != vReplaceCells.end(); ++c_iter) {
+          const ALE::Obj<sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, depth);
+
+          if (preFace->size() == faceSize) {
+            if (debug) {std::cout << "    Scheduling " << *n_iter << " for replacement" << std::endl;}
+            vReplaceCells.insert(*n_iter);
+            modified   = true;
+            classified = true;
+            break;
+          }
+        }
+        if (classified) continue;
+        // It is unclear whether taking out the noReplace cells will speed this up
+        for(PointSet::const_iterator c_iter = vNoReplaceCells.begin(); c_iter != vNoReplaceCells.end(); ++c_iter) {
+          const ALE::Obj<sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, depth);
+
+          if (preFace->size() == faceSize) {
+            if (debug) {std::cout << "    Scheduling " << *n_iter << " for no replacement" << std::endl;}
+            vNoReplaceCells.insert(*n_iter);
+            modified   = true;
+            classified = true;
+            break;
+          }
+        }
+      }
+      if (debug) {
+        std::cout << "classifySize: " << classifySize << std::endl;
+        std::cout << "classifyTotal: " << classifyTotal << std::endl;
+        std::cout << "vReplaceCells.size: " << vReplaceCells.size() << std::endl;
+        std::cout << "vNoReplaceCells.size: " << vNoReplaceCells.size() << std::endl;
+      }
+      assert(classifySize < vReplaceCells.size() + vNoReplaceCells.size());
+      classifySize = vReplaceCells.size() + vNoReplaceCells.size();
+      assert(classifySize <= classifyTotal);
+    }
+    replaceCells.insert(vReplaceCells.begin(), vReplaceCells.end());
+  }
+  debug = 0;
+#else
+  for(PointSet::const_iterator v_iter = replaceVertices.begin(); v_iter != replaceVertices.end(); ++v_iter) {
     bool modified = true;
 
+    // Replace cells that share a face with this one, or with one that does
     while(modified) {
       modified = false;
       const ALE::Obj<sieve_type::traits::supportSequence>& neighbors = sieve->support(*v_iter);
@@ -431,7 +513,7 @@
         if (debug) std::cout << "  Checking fault neighbor " << *n_iter << std::endl;
         // If neighbors shares a faces with anyone in replaceCells, then add
         for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
-          const ALE::Obj<sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, mesh->depth());
+          const ALE::Obj<sieve_type::coneSet>& preFace = sieve->nMeet(*c_iter, *n_iter, depth);
 
           if (preFace->size() == faceSize) {
             if (debug) std::cout << "    Scheduling " << *n_iter << " for replacement" << std::endl;
@@ -443,6 +525,7 @@
       }
     }
   }
+#endif
   for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
     _replaceCell(sieve, *c_iter, &vertexRenumber, debug);
   }

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc	2007-09-26 20:12:18 UTC (rev 8034)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc	2007-09-27 03:32:37 UTC (rev 8035)
@@ -575,7 +575,7 @@
   ALE::Obj<ALE::Mesh> mesh;
   meshio::MeshIOAscii iohandler;
   iohandler.filename(data.filename);
-  iohandler.debug(true);
+  iohandler.debug(false);
   iohandler.interpolate(false);
   iohandler.read(&mesh);
 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataQuad4h.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataQuad4h.cc	2007-09-26 20:12:18 UTC (rev 8034)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataQuad4h.cc	2007-09-27 03:32:37 UTC (rev 8035)
@@ -37,21 +37,21 @@
  *
  * 9 ----10 ----25-11 ----12
  * |      |      |  |      |
+ * |  0   |  1   |28|  2   |
  * |      |      |  |      |
  * |      |      |  |      |
- * |      |      |  |      |
  *13 ----14 ----26-15 ----16
  * |      |      |  |      |
+ * |  3   |  4   |29|  5   |
  * |      |      |  |      |
  * |      |      |  |      |
- * |      |      |  |      |
- *17 ----18 ----27-19 ----20
+ *30 ----31 ----27-19 ----20
+ * |  32  |      |         |
+ *17 ----18      |         |
+ * |      |   7  |    8    |
+ * |  6   |      |         |
  * |      |      |         |
- *28 ----29      |         |
  * |      |      |         |
- * |      |      |         |
- * |      |      |         |
- * |      |      |         |
  *21 ----22 ----23 -------24
  *
  */
@@ -113,8 +113,8 @@
   14, 31, 27, 26,
   15, 19, 20, 16,
   17, 21, 22, 18,
-  31, 22, 23, 27,
-  27, 23, 24, 20,
+  18, 22, 23, 27,
+  19, 23, 24, 20,
   11, 15, 25, 26,
   15, 19, 26, 27,
   18, 17, 31, 30,

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataTri3f.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataTri3f.cc	2007-09-26 20:12:18 UTC (rev 8034)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataTri3f.cc	2007-09-27 03:32:37 UTC (rev 8035)
@@ -15,14 +15,14 @@
  * Cells are 0-3, vertices are 4-9.
  *
  *       5 -------- 7 -------- 9
- *       |        / | \        |
- *       |       /  |  \       |
- *       |      /   |   \      |
- *       |     /    |    \     |
- *       |    /     |     \    |
- *       |   /      |      \   |
- *       |  /       |       \  |
- *       | /        |        \ |
+ *       |        / |        / |
+ *       |       /  |       /  |
+ *       |  0   /   |  1   /   |
+ *       |     /    |     /    |
+ *       |    /     |    /     |
+ *       |   /      |   /      |
+ *       |  /   3   |  /   2   |
+ *       | /        | /        |
  *       4 -------- 6 -------- 8
  *
  * After adding cohesive elements
@@ -32,13 +32,13 @@
  *       5 -------- 7 --11 -------- 9
  *       |        / |    |        / |
  *       |       /  |    |       /  |
- *       |      /   |    |      /   |
- *       |     /    |    |     /    |
+ *       |  0   /   |    |  1   /   |
+ *       |     /    | 12 |     /    |
  *       |    /     |    |    /     |
  *       |   /      |    |   /      |
- *       |  /       |    |  /       |
+ *       |  /   3   |    |  /   2   |
  *       | /        |    | /        |
- *       4 -------- 6 -- 6 -------- 8
+ *       4 -------- 6 --10 -------- 8
  */
 
 #include "CohesiveDataTri3f.hh"



More information about the cig-commits mailing list