diff --git a/dune/faultnetworks/faultfactories/rockfaultfactory.hh b/dune/faultnetworks/faultfactories/rockfaultfactory.hh index e8f216e60b7295bc40f132fe125c9f5a7d3efcad..432cc445e8cf50f851ec708d0c958c48eff71e0a 100644 --- a/dune/faultnetworks/faultfactories/rockfaultfactory.hh +++ b/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); } }