.

dune/faultnetworks/faultfactories/rockfaultfactory.hh

@@ -29,15 +29,9 @@ public:
     int top;
     int bottom;
 
-    // fault idx
-    int leftFault;
-    int rightFault;
-    int topFault;
-    int bottomFault;
-
     Dune::FieldVector<ctype, 2> center;
 
-    void setVertices(int _left, int _right, int _top, int _bottom, const Dune::FieldVector<ctype, 2>& _center) {
+    void set(int _left, int _right, int _top, int _bottom, const Dune::FieldVector<ctype, 2>& _center) {
         left = _left;
         right = _right;
         top = _top;
@@ -45,13 +39,6 @@ public:
 
         center = _center;
     }
-
-    void setFaults(int _left, int _right, int _top, int _bottom) {
-        leftFault = _left;
-        rightFault = _right;
-        topFault = _top;
-        bottomFault = _bottom;
-    }
 };
 
 
@@ -442,6 +429,18 @@ private:
         return lastDof;
     }
 
+    void createRock(Rock& 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) {
+
+        rock.level = level_;
+        rock.center = center;
+
+        const ID centerIDs = {computeID(center(), 0), computeID(center(), 1)};
+
+        rock.set(left, right, top, bottom, center);
+    }
+
     void prolong(const Rock& rock, Rock& newRock) {
         newRock.level = rock.level;
 
@@ -450,8 +449,7 @@ private:
         auto newTop = coarseToLevelVertex_[rock.top()];
         auto newBottom = coarseToLevelVertex_[rock.bottom()];
 
-        newRock.setVertices(newLeft, newRight, newTop, newBottom, rock.center());
-        newRock.setFaults(rock.leftFault, rock.rightFault, rock.topFault, rock.bottomFault);
+        newRock.set(newLeft, newRight, newTop, newBottom, rock.center());
     }
 
     bool randomSplit(const Rock& rock) {
@@ -479,29 +477,54 @@ private:
         if (!toBeSplit) {
             rocks_.push_back(newRock);
         } else {
-            const ID centerIDs = {computeID(newRock.center(), 0), computeID(newRock.center(), 1)};
+            const ID centerIDs = {computeID(newRock.center, 0), computeID(newRock.center, 1)};
 
-            size_t xCorridor = 1.0/2 * std::min(centerIDs[0] - vertexIDs_[newRock.left()][0],
-                    vertexIDs_[newRock.right()][0] - centerIDs[0]) + 1;
-            size_t yCorridor = 1.0/2 * std::min(centerIDs[1] - vertexIDs_[newRock.bottom()][1],
-                    vertexIDs_[newRock.top()][1] - centerIDs[1]) + 1;
+            size_t xCorridor = 1.0/2 * std::min(centerIDs[0] - vertexIDs_[newRock.left][0],
+                    vertexIDs_[newRock.right][0] - centerIDs[0]) + 1;
+            size_t yCorridor = 1.0/2 * std::min(centerIDs[1] - vertexIDs_[newRock.bottom][1],
+                    vertexIDs_[newRock.top][1] - centerIDs[1]) + 1;
 
             // split rock into 4 subparts by 4 new faults intersecting at center of rock
             std::shared_ptr<FaultInterface<GV>> x1Fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
-            generateFault(x1Fault, newRock.left(), newRock.center(), yCorridor, separatingDofs);
+            generateFault(x1Fault, newRock.left, newRock.center, yCorridor, separatingDofs);
             faults_.push_back(x1Fault);
 
             std::shared_ptr<FaultInterface<GV>> x2Fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
-            generateFault(x2Fault, newRock.right(), newRock.center(), yCorridor, separatingDofs);
+            generateFault(x2Fault, newRock.right, newRock.center, yCorridor, separatingDofs);
             faults_.push_back(x2Fault);
 
             std::shared_ptr<FaultInterface<GV>> y1Fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
-            generateFault(y1Fault, newRock.bottom(), newRock.center(), xCorridor, separatingDofs);
+            generateFault(y1Fault, newRock.bottom, newRock.center, xCorridor, separatingDofs);
             faults_.push_back(y1Fault);
 
             std::shared_ptr<FaultInterface<GV>> y2Fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
-            generateFault(y2Fault, newRock.top(), newRock.center(), xCorridor, separatingDofs);
+            generateFault(y2Fault, newRock.top, newRock.center, xCorridor, separatingDofs);
             faults_.push_back(y2Fault);
+
+            auto left = vertexPositions_[newRock.left];
+            auto right = vertexPositions_[newRock.right];
+            auto top = vertexPositions_[newRock.top];
+            auto bottom = vertexPositions_[newRock.bottom];
+
+            Rock rock00;
+            auto center00 = 1.0/2*(right + bottom);
+            createRock(rock00, center00, separatingDofs, x2Fault.getInterfaceDofs(), 0, y1Fault.getInterfaceDofs(), 0);
+            rocks_.push_back(rock00);
+
+            Rock rock01;
+            auto center01 = 1.0/2*(left + bottom);
+            createRock(rock01, center01, separatingDofs, x1Fault.getInterfaceDofs(), 0, y1Fault.getInterfaceDofs(), 1);
+            rocks_.push_back(rock01);
+
+            Rock rock10;
+            auto center10 = 1.0/2*(right + top);
+            createRock(rock10, center10, separatingDofs, x2Fault.getInterfaceDofs(), 1, y2Fault.getInterfaceDofs(), 0);
+            rocks_.push_back(rock10);
+
+            Rock rock11;
+            auto center11 = 1.0/2*(left + top);
+            createRock(rock11, center11, separatingDofs, x1Fault.getInterfaceDofs(), 1, y2Fault.getInterfaceDofs(), 1);
+            rocks_.push_back(rock11);
         }
     }