Skip to content
Snippets Groups Projects
Commit 8cc655c0 authored by podlesny's avatar podlesny
Browse files

.

parent 2ab1e070
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/faultnetworks/faultfactories/rockfaultfactory.hh
+ 127
− 49
View file @ 8cc655c0
...@@ -20,50 +20,37 @@ ...@@ -20,50 +20,37 @@
#include <dune/fufem/boundaryiterator.hh> #include <dune/fufem/boundaryiterator.hh>
template<class ctype = double> class Rock { 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: public:
void set(int level, int left, int right, int top, int bottom, const Dune::FieldVector<ctype, 2>& center) { int level;
level_ = level;
left_ = left;
right_ = right;
top_ = top;
bottom_ = bottom;
center_ = center;
}
auto level() const { // vertex idx
return level_; int left;
} int right;
int top;
int bottom;
auto left() const { // fault idx
return left_; int leftFault;
} int rightFault;
int topFault;
int bottomFault;
auto right() const { Dune::FieldVector<ctype, 2> center;
return right_;
}
auto top() const { void setVertices(int _left, int _right, int _top, int _bottom, const Dune::FieldVector<ctype, 2>& _center) {
return top_; left = _left;
} right = _right;
top = _top;
bottom = _bottom;
auto bottom() const { center = _center;
return bottom_;
} }
const auto& center() { void setFaults(int _left, int _right, int _top, int _bottom) {
return center_; leftFault = _left;
rightFault = _right;
topFault = _top;
bottomFault = _bottom;
} }
}; };
...@@ -102,6 +89,8 @@ protected: ...@@ -102,6 +89,8 @@ protected:
const int level_; const int level_;
const ctype resolution_; const ctype resolution_;
const GV gridView_; const GV gridView_;
const double splittingThreshold_;
const double maxAngle_; const double maxAngle_;
const typename GV::IndexSet& indexSet_; const typename GV::IndexSet& indexSet_;
...@@ -110,6 +99,11 @@ protected: ...@@ -110,6 +99,11 @@ protected:
std::vector<Intersection> faces_; std::vector<Intersection> faces_;
std::vector<Coords> vertexPositions_; 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<std::vector<size_t>> vertexToFaces_;
std::vector<int> coarseToLevelVertex_; std::vector<int> coarseToLevelVertex_;
...@@ -285,11 +279,11 @@ private: ...@@ -285,11 +279,11 @@ private:
dim = 1; 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; std::set<size_t> separatingIDs;
auto faultSeedID = computeID(vertexPositions_[faultSeedIdx], dim); auto faultSeedID = vertexIDs_[faultSeedIdx][dim];
separatingIDs.insert(faultSeedID + corridor); separatingIDs.insert(faultSeedID + corridor);
separatingIDs.insert(faultSeedID - corridor); separatingIDs.insert(faultSeedID - corridor);
...@@ -380,13 +374,84 @@ private: ...@@ -380,13 +374,84 @@ private:
//std::cout << "------------------------------------- " << std::endl << std::endl; //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) { void prolong(const Rock& rock, Rock& newRock) {
newRock.level = rock.level;
auto newLeft = coarseToLevelVertex_[rock.left()]; auto newLeft = coarseToLevelVertex_[rock.left()];
auto newRight = coarseToLevelVertex_[rock.right()]; auto newRight = coarseToLevelVertex_[rock.right()];
auto newTop = coarseToLevelVertex_[rock.top()]; auto newTop = coarseToLevelVertex_[rock.top()];
auto newBottom = coarseToLevelVertex_[rock.bottom()]; 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) { bool randomSplit(const Rock& rock) {
...@@ -394,29 +459,32 @@ private: ...@@ -394,29 +459,32 @@ private:
double res = 0.0; double res = 0.0;
if (center[0] < center[1]) { if (center[0] < center[1]) {
res = 100.0*(1.0 - center[0]); res = 1.0 - center[0];
} else { } 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; 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) { void split(const Rock& rock, const std::set<size_t>& separatingDofs) {
Rock newRock; Rock newRock;
prolong(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) { if (!toBeSplit) {
rocks_.push_back(newRock); rocks_.push_back(newRock);
} else { } 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), size_t xCorridor = 1.0/2 * std::min(centerIDs[0] - vertexIDs_[newRock.left()][0],
computeID(vertexPositions_[newRock.right()], 0) - centerIDs[0]) + 1; vertexIDs_[newRock.right()][0] - centerIDs[0]) + 1;
size_t yCorridor = 1.0/2 * std::min(centerIDs[1] - computeID(vertexPositions_[newRock.bottom()], 1), size_t yCorridor = 1.0/2 * std::min(centerIDs[1] - vertexIDs_[newRock.bottom()][1],
computeID(vertexPositions_[newRock.top()], 1) - centerIDs[1]) + 1; vertexIDs_[newRock.top()][1] - centerIDs[1]) + 1;
// split rock into 4 subparts by 4 new faults intersecting at center of rock // 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_); std::shared_ptr<FaultInterface<GV>> x1Fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
...@@ -450,6 +518,7 @@ private: ...@@ -450,6 +518,7 @@ private:
faces_.resize(faceMapper_.size()); faces_.resize(faceMapper_.size());
vertexPositions_.resize(gridView_.size(dim)); vertexPositions_.resize(gridView_.size(dim));
vertexToFaces_.resize(gridView_.size(dim)); vertexToFaces_.resize(gridView_.size(dim));
vertexIDs_.resize(gridView_.size(dim));
coarseToLevelVertex_.resize(coarseLevelFactory_.gridView().size(dim)); coarseToLevelVertex_.resize(coarseLevelFactory_.gridView().size(dim));
std::vector<bool> faceHandled(faceMapper_.size(), false); std::vector<bool> faceHandled(faceMapper_.size(), false);
...@@ -476,9 +545,14 @@ private: ...@@ -476,9 +545,14 @@ private:
size_t idxInElement = refElement.subEntity(isect.indexInInside(), 1, i, dimElement); size_t idxInElement = refElement.subEntity(isect.indexInInside(), 1, i, dimElement);
const size_t vertexID = indexSet_.subIndex(elem, idxInElement, 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); vertexToFaces_[vertexID].push_back(faceID);
ID id = {computeID(vertex, 0), computeID(vertex, 1)};
vertexIDs_[vertexID] = id;
IDsToDof_[id] = vertexID;
} }
} }
} }
...@@ -509,6 +583,10 @@ private: ...@@ -509,6 +583,10 @@ private:
const auto& rocks() const { const auto& rocks() const {
return rocks_; return rocks_;
} }
auto level() const {
return level_;
}
}; };
template <class GridType> template <class GridType>
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment