.

dune/faultnetworks/faultfactories/rockfaultfactory.hh

@@ -70,7 +70,7 @@ protected:
     using FaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout >;
 
 
-    using Rock = Rock<ctype>;
+    using MyRock = Rock<ctype>;
 
     const std::shared_ptr<GridType> grid_;
     const int level_;
@@ -95,7 +95,7 @@ protected:
     std::vector<int> coarseToLevelVertex_;
 
     const LevelRockFaultFactory& coarseLevelFactory_;
-    std::vector<Rock> rocks_;
+    std::vector<MyRock> rocks_;
     std::vector<std::shared_ptr<FaultInterface<GV>>>& faults_;
 
 private:
@@ -363,73 +363,30 @@ private:
 
     auto searchDof(const ID& IDs, const std::set<size_t>& separatingDofs, size_t dim, int dir) {
         auto candidatIDs = IDs;
-        size_t count = 0;
-
         int lastDof = -1;
-        bool keepSearching = true;
 
         if (dir>0) {
             while (lastDof<0) {
-                candidatIDs[dim] = IDs[dim] - count;
+                candidatIDs[dim]++;
                 auto dof = IDsToDof_[candidatIDs];
                 if (separatingDofs.count(dof)) {
                     lastDof = dof;
                 }
-
-                candidatIDs[dim] = IDs[dim] + count;
-                dof = IDsToDof_[candidatIDs];
-                if (separatingDofs.count(dof)) {
-                    lastDof = dof;
-                }
-
-                count++;
-            }
-
-            while (keepSearching) {
-                candidatIDs[dim] = IDs[dim] + count;
-                auto dof = IDsToDof_[candidatIDs];
-                if (separatingDofs.count(dof)) {
-                    lastDof = dof;
-                } else {
-                    keepSearching = false;
-                }
-
-                count ++;
             }
         } else {
             while (lastDof<0) {
-                candidatIDs[dim] = IDs[dim] + count;
-                auto dof = IDsToDof_[candidatIDs];
-                if (separatingDofs.count(dof)) {
-                    lastDof = dof;
-                }
-
-                candidatIDs[dim] = IDs[dim] - count;
-                dof = IDsToDof_[candidatIDs];
-                if (separatingDofs.count(dof)) {
-                    lastDof = dof;
-                }
-
-                count++;
-            }
-
-            while (keepSearching) {
-                candidatIDs[dim] = IDs[dim] - count;
+                candidatIDs[dim]--;
                 auto dof = IDsToDof_[candidatIDs];
                 if (separatingDofs.count(dof)) {
                     lastDof = dof;
-                } else {
-                    keepSearching = false;
                 }
-
-                count ++;
             }
         }
 
         return lastDof;
     }
 
-    void createRock(Rock& rock, const Coords& center, const std::set<size_t>& separatingDofs,
+    void createRock(MyRock& rock, const Coords& center, const std::set<size_t>& separatingDofs,
                     const std::set<size_t>& xFaultDofs, bool xFaultBottom,
                     const std::set<size_t>& yFaultDofs, bool yFaultRight) {
 
@@ -438,10 +395,28 @@ private:
 
         const ID centerIDs = {computeID(center(), 0), computeID(center(), 1)};
 
+        size_t left, right, top, bottom = 0;
+
+        if (xFaultBottom) {
+            top = searchDof(centerIDs, separatingDofs, 1, 1);
+            bottom = searchDof(centerIDs, xFaultDofs, 1, -1);
+        } else {
+            top = searchDof(centerIDs, xFaultDofs, 1, 1);
+            bottom = searchDof(centerIDs, separatingDofs, 1, -1);
+        }
+
+        if (yFaultRight) {
+            left = searchDof(centerIDs, separatingDofs, 0, -1);
+            right = searchDof(centerIDs, yFaultDofs, 0, 1);
+        } else {
+            left = searchDof(centerIDs, yFaultDofs, 0, -1);
+            right = searchDof(centerIDs, separatingDofs, 0, 1);
+        }
+
         rock.set(left, right, top, bottom, center);
     }
 
-    void prolong(const Rock& rock, Rock& newRock) {
+    void prolong(const MyRock& rock, MyRock& newRock) {
         newRock.level = rock.level;
 
         auto newLeft = coarseToLevelVertex_[rock.left()];
@@ -452,7 +427,7 @@ private:
         newRock.set(newLeft, newRight, newTop, newBottom, rock.center());
     }
 
-    bool randomSplit(const Rock& rock) {
+    bool randomSplit(const MyRock& rock) {
         const auto& center = rock.center();
 
         double res = 0.0;
@@ -469,7 +444,7 @@ private:
         return  100.0*res > prob;
     }
 
-    void split(const Rock& rock, const std::set<size_t>& separatingDofs) {
+    void split(const MyRock& rock, const std::set<size_t>& separatingDofs) {
         Rock newRock;
         prolong(rock, newRock);
 
@@ -506,22 +481,22 @@ private:
             auto top = vertexPositions_[newRock.top];
             auto bottom = vertexPositions_[newRock.bottom];
 
-            Rock rock00;
+            MyRock rock00;
             auto center00 = 1.0/2*(right + bottom);
             createRock(rock00, center00, separatingDofs, x2Fault.getInterfaceDofs(), 0, y1Fault.getInterfaceDofs(), 0);
             rocks_.push_back(rock00);
 
-            Rock rock01;
+            MyRock rock01;
             auto center01 = 1.0/2*(left + bottom);
             createRock(rock01, center01, separatingDofs, x1Fault.getInterfaceDofs(), 0, y1Fault.getInterfaceDofs(), 1);
             rocks_.push_back(rock01);
 
-            Rock rock10;
+            MyRock rock10;
             auto center10 = 1.0/2*(right + top);
             createRock(rock10, center10, separatingDofs, x2Fault.getInterfaceDofs(), 1, y2Fault.getInterfaceDofs(), 0);
             rocks_.push_back(rock10);
 
-            Rock rock11;
+            MyRock rock11;
             auto center11 = 1.0/2*(left + top);
             createRock(rock11, center11, separatingDofs, x1Fault.getInterfaceDofs(), 1, y2Fault.getInterfaceDofs(), 1);
             rocks_.push_back(rock11);