diff --git a/dune/faultnetworks/faultfactories/rockfaultfactory.hh b/dune/faultnetworks/faultfactories/rockfaultfactory.hh
index 432cc445e8cf50f851ec708d0c958c48eff71e0a..b240b7a2af9ce20029d897a29f6e89c8fd161104 100644
--- a/dune/faultnetworks/faultfactories/rockfaultfactory.hh
+++ b/dune/faultnetworks/faultfactories/rockfaultfactory.hh
@@ -20,50 +20,37 @@
#include <dune/fufem/boundaryiterator.hh>
template<class ctype = double> class Rock {
-protected:
- int level_;
-
- // vertex IDs
- int left_;
- int right_;
- int top_;
- int bottom_;
-
- Dune::FieldVector<ctype, 2> center_;
public:
- void set(int level, int left, int right, int top, int bottom, const Dune::FieldVector<ctype, 2>& center) {
- level_ = level;
-
- left_ = left;
- right_ = right;
- top_ = top;
- bottom_ = bottom;
-
- center_ = center;
- }
+ int level;
- auto level() const {
- return level_;
- }
+ // vertex idx
+ int left;
+ int right;
+ int top;
+ int bottom;
- auto left() const {
- return left_;
- }
+ // fault idx
+ int leftFault;
+ int rightFault;
+ int topFault;
+ int bottomFault;
- auto right() const {
- return right_;
- }
+ Dune::FieldVector<ctype, 2> center;
- auto top() const {
- return top_;
- }
+ void setVertices(int _left, int _right, int _top, int _bottom, const Dune::FieldVector<ctype, 2>& _center) {
+ left = _left;
+ right = _right;
+ top = _top;
+ bottom = _bottom;
- auto bottom() const {
- return bottom_;
+ center = _center;
}
- const auto& center() {
- return center_;
+ void setFaults(int _left, int _right, int _top, int _bottom) {
+ leftFault = _left;
+ rightFault = _right;
+ topFault = _top;
+ bottomFault = _bottom;
}
};
@@ -102,6 +89,8 @@ protected:
const int level_;
const ctype resolution_;
const GV gridView_;
+
+ const double splittingThreshold_;
const double maxAngle_;
const typename GV::IndexSet& indexSet_;
@@ -110,6 +99,11 @@ protected:
std::vector<Intersection> faces_;
std::vector<Coords> vertexPositions_;
+
+ using ID = std::array<size_t, 2>;
+ std::vector<ID> vertexIDs_;
+ std::map<ID, size_t> IDsToDof_;
+
std::vector<std::vector<size_t>> vertexToFaces_;
std::vector<int> coarseToLevelVertex_;
@@ -285,11 +279,11 @@ private:
dim = 1;
}
- const std::array<size_t, 2> centerIDs = {computeID(center, 0), computeID(center, 1)};
+ const ID centerIDs = {computeID(center, 0), computeID(center, 1)};
std::set<size_t> separatingIDs;
- auto faultSeedID = computeID(vertexPositions_[faultSeedIdx], dim);
+ auto faultSeedID = vertexIDs_[faultSeedIdx][dim];
separatingIDs.insert(faultSeedID + corridor);
separatingIDs.insert(faultSeedID - corridor);
@@ -380,13 +374,84 @@ private:
//std::cout << "------------------------------------- " << std::endl << std::endl;
}
+ 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;
+ 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;
+ auto dof = IDsToDof_[candidatIDs];
+ if (separatingDofs.count(dof)) {
+ lastDof = dof;
+ } else {
+ keepSearching = false;
+ }
+
+ count ++;
+ }
+ }
+
+ return lastDof;
+ }
+
void prolong(const Rock& rock, Rock& newRock) {
+ newRock.level = rock.level;
+
auto newLeft = coarseToLevelVertex_[rock.left()];
auto newRight = coarseToLevelVertex_[rock.right()];
auto newTop = coarseToLevelVertex_[rock.top()];
auto newBottom = coarseToLevelVertex_[rock.bottom()];
- newRock.set(rock.level(), newLeft, newRight, newTop, newBottom, rock.center());
+ newRock.setVertices(newLeft, newRight, newTop, newBottom, rock.center());
+ newRock.setFaults(rock.leftFault, rock.rightFault, rock.topFault, rock.bottomFault);
}
bool randomSplit(const Rock& rock) {
@@ -394,29 +459,32 @@ private:
double res = 0.0;
if (center[0] < center[1]) {
- res = 100.0*(1.0 - center[0]);
+ res = 1.0 - center[0];
} else {
- res = 100.0*(1.0 - center[1]);
+ res = 1.0 - center[1];
}
+ if (res > splittingThreshold_)
+ return true;
+
double prob = (100.0 * std::rand() / (RAND_MAX + 1.0)) + 1;
- return res > prob;
+ return 100.0*res > prob;
}
void split(const Rock& rock, const std::set<size_t>& separatingDofs) {
Rock newRock;
prolong(rock, newRock);
- bool toBeSplit = (rock.level() == level_-1) and randomSplit(rock);
+ bool toBeSplit = (rock.level() == coarseLevelFactory_.level()) and randomSplit(rock);
if (!toBeSplit) {
rocks_.push_back(newRock);
} else {
- const std::array<size_t, 2> 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] - 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;
+ 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_);
@@ -450,6 +518,7 @@ private:
faces_.resize(faceMapper_.size());
vertexPositions_.resize(gridView_.size(dim));
vertexToFaces_.resize(gridView_.size(dim));
+ vertexIDs_.resize(gridView_.size(dim));
coarseToLevelVertex_.resize(coarseLevelFactory_.gridView().size(dim));
std::vector<bool> faceHandled(faceMapper_.size(), false);
@@ -476,9 +545,14 @@ private:
size_t idxInElement = refElement.subEntity(isect.indexInInside(), 1, i, dimElement);
const size_t vertexID = indexSet_.subIndex(elem, idxInElement, dimElement);
+ const auto& vertex = elem.geometry().corner(idxInElement);
- vertexPositions_[vertexID] = elem.geometry().corner(idxInElement);
+ vertexPositions_[vertexID] = vertex;
vertexToFaces_[vertexID].push_back(faceID);
+
+ ID id = {computeID(vertex, 0), computeID(vertex, 1)};
+ vertexIDs_[vertexID] = id;
+ IDsToDof_[id] = vertexID;
}
}
}
@@ -509,6 +583,10 @@ private:
const auto& rocks() const {
return rocks_;
}
+
+ auto level() const {
+ return level_;
+ }
};
template <class GridType>