.

dune/faultnetworks/faultfactories/rockfaultfactory.hh

@@ -124,6 +124,7 @@ private:
 
         bool intersectionAllowed(const size_t faceID, const size_t vertexID,
                                  const std::array<size_t, 2>& centerIDs,
+                                 const std::set<size_t>& deadendDofs,
                                  const std::set<int>& separatingIDs, const std::set<size_t>& faultDofs,
                                  const Coords& direction, size_t dim) const {
 
@@ -148,7 +149,7 @@ private:
                 std::array<size_t, 2> IDs = {computeID(vertex, 0), computeID(vertex, 1)};
 
                 bool centerPassed = (sgn(direction[dim])>0) ? (centerIDs[dim]<IDs[dim]) : (centerIDs[dim]>IDs[dim]);
-                if (faultDofs.count(isDof) or separatingIDs.count(IDs[otherDim]) or centerPassed) {
+                if (faultDofs.count(isDof) or separatingIDs.count(IDs[otherDim]) or centerPassed or deadendDofs.count(isDof)) {
                     return false;
                 }
             }
@@ -267,13 +268,13 @@ private:
 
         void generateFault(std::shared_ptr<FaultInterface<GV>> fault, const size_t faultSeedIdx, const Coords& center,
                            const size_t corridor,
-                           const std::set<size_t>& constSeparatingDofs) {
+                           const std::set<size_t>& separatingDofs) {
             //std::cout << "LevelGeoFaultFactory::generateFault() " << std::endl;
             //std::cout << "------------------------------------- " << std::endl << std::endl;
 
             bool success = false;
 
-            std::set<size_t> separatingDofs = constSeparatingDofs;
+            std::set<size_t> deadendDofs = separatingDofs;
 
             auto direction = center;
             direction -= vertexPositions_[faultSeedIdx];
@@ -316,7 +317,7 @@ private:
                     std::vector<size_t> admissibleFaces;
 
                     for (size_t i=0; i<faces.size(); i++) {
-                        if (intersectionAllowed(faces[i], vertexID, centerIDs, separatingDofs, separatingIDs, faultDofs, direction, dim)) {
+                        if (intersectionAllowed(faces[i], vertexID, centerIDs, deadendDofs, separatingIDs, faultDofs, direction, dim)) {
                             admissibleFaces.push_back(faces[i]);
                         }
                     }
@@ -339,7 +340,7 @@ private:
                     faultFaces.pop_back();
                     vertexToAdmissibleFaces.erase(vertexID);
 
-                    separatingDofs.insert(vertexID);
+                    deadendDofs.insert(vertexID);
                     vertexQueue.push(faultDofs.back());
                 } else {
 
@@ -410,14 +411,29 @@ private:
         if (!toBeSplit) {
             rocks_.push_back(newRock);
         } else {
-            // split rock into 4 subparts by two new intersecting faults
-            std::shared_ptr<FaultInterface<GV>> xFault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
-            generateFault(xFault, newRock, separatingDofs, 0);
-            faults_.push_back(xFault);
-
-            std::shared_ptr<FaultInterface<GV>> yFault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
-            generateFault(yFault, newRock, separatingDofs, 1);
-            faults_.push_back(yFault);
+            const std::array<size_t, 2> centerIDs = {computeID(newRock.center(), 0), computeID(newRock.center(), 1)};
+
+            size_t xCorridor = 1.0/2 * std::min(centerIDs[0] - computeID(vertexPositions_[newRock.left()], 0),
+                    computeID(vertexPositions_[newRock.right()], 0) - centerIDs[0]) + 1;
+            size_t yCorridor = 1.0/2 * std::min(centerIDs[1] - computeID(vertexPositions_[newRock.bottom()], 1),
+                    computeID(vertexPositions_[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);
+            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);
+            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);
+            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);
+            faults_.push_back(y2Fault);
         }
     }