diff --git a/dune/faultnetworks/faultfactories/rockfaultfactory.hh b/dune/faultnetworks/faultfactories/rockfaultfactory.hh index 82d8b0d651ea89f71d641f2eef1a81994cdaaecc..e48d56369d744a0efb4fd44527384ae55f08a8f2 100644 --- a/dune/faultnetworks/faultfactories/rockfaultfactory.hh +++ b/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);