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

.

parent b33d8c8f
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
dune/faultnetworks/faultfactories/rockfaultfactory.hh
+ 67
− 51
View file @ b0c4be5c
...@@ -91,22 +91,20 @@ protected: ...@@ -91,22 +91,20 @@ protected:
std::vector<ID> vertexIDs_; std::vector<ID> vertexIDs_;
std::map<ID, size_t> IDsToDof_; std::map<ID, size_t> IDsToDof_;
std::set<size_t> boundaryDofs_;
std::vector<std::vector<size_t>> vertexToFaces_; std::vector<std::vector<size_t>> vertexToFaces_;
std::vector<int> coarseToLevelVertex_; std::vector<int> coarseToLevelVertex_;
const LevelRockFaultFactory& coarseLevelFactory_; const LevelRockFaultFactory& coarseLevelFactory_;
std::vector<MyRock> rocks_; std::vector<MyRock> rocks_;
std::vector<std::shared_ptr<FaultInterface<GV>>>& faults_; std::vector<std::shared_ptr<FaultInterface<GV>>> faults_;
private:
template <typename T> int sgn(T val) {
return (T(0) < val) - (val < T(0));
}
protected:
bool intersectionAllowed(const size_t faceID, const size_t vertexID, bool intersectionAllowed(const size_t faceID, const size_t vertexID,
const std::array<size_t, 2>& centerIDs, const std::array<size_t, 2>& centerIDs,
const std::set<size_t>& deadendDofs, const std::set<size_t>& deadendDofs,
const std::set<int>& separatingIDs, const std::set<size_t>& faultDofs, const std::set<size_t>& separatingIDs, const std::set<size_t>& faultDofs,
const Coords& direction, size_t dim) const { const Coords& direction, size_t dim) const {
const auto& intersection = faces_[faceID]; const auto& intersection = faces_[faceID];
...@@ -127,9 +125,10 @@ private: ...@@ -127,9 +125,10 @@ private:
const size_t otherDim = (dim + 1) % 2; const size_t otherDim = (dim + 1) % 2;
for (const auto& isDof : intersectionDofs) { for (const auto& isDof : intersectionDofs) {
const auto& vertex = vertexPositions_[isDof]; const auto& vertex = vertexPositions_[isDof];
std::array<size_t, 2> IDs = {computeID(vertex, 0), computeID(vertex, 1)}; ID IDs = {computeID(vertex, 0), computeID(vertex, 1)};
bool centerPassed = (sgn(direction[dim])>0) ? (centerIDs[dim]<IDs[dim]) : (centerIDs[dim]>IDs[dim]); auto val = direction[dim];
bool centerPassed = ((0 < val) - (val < 0)>0) ? (centerIDs[dim]<IDs[dim]) : (centerIDs[dim]>IDs[dim]);
if (faultDofs.count(isDof) or separatingIDs.count(IDs[otherDim]) or centerPassed or deadendDofs.count(isDof)) { if (faultDofs.count(isDof) or separatingIDs.count(IDs[otherDim]) or centerPassed or deadendDofs.count(isDof)) {
return false; return false;
} }
...@@ -250,17 +249,19 @@ private: ...@@ -250,17 +249,19 @@ private:
if (std::abs(direction[0]) < std::abs(direction[1])) { if (std::abs(direction[0]) < std::abs(direction[1])) {
dim = 1; dim = 1;
} }
auto otherDim = (dim+1) % 2;
const ID 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 = vertexIDs_[faultSeedIdx][dim]; auto faultSeedID = vertexIDs_[faultSeedIdx][otherDim];
separatingIDs.insert(faultSeedID + corridor); separatingIDs.insert(faultSeedID + corridor);
separatingIDs.insert(faultSeedID - corridor); separatingIDs.insert(faultSeedID - corridor);
std::vector<size_t> faultDofs(0); std::vector<size_t> faultDofs(0);
std::set<size_t> faultDofSet;
faultDofs.push_back(faultSeedIdx); faultDofs.push_back(faultSeedIdx);
faultDofSet.insert(faultSeedIdx);
std::vector<size_t> faultFaces; std::vector<size_t> faultFaces;
...@@ -283,7 +284,7 @@ private: ...@@ -283,7 +284,7 @@ private:
std::vector<size_t> admissibleFaces; std::vector<size_t> admissibleFaces;
for (size_t i=0; i<faces.size(); i++) { for (size_t i=0; i<faces.size(); i++) {
if (intersectionAllowed(faces[i], vertexID, centerIDs, deadendDofs, separatingIDs, faultDofs, direction, dim)) { if (intersectionAllowed(faces[i], vertexID, centerIDs, deadendDofs, separatingIDs, faultDofSet, direction, dim)) {
admissibleFaces.push_back(faces[i]); admissibleFaces.push_back(faces[i]);
} }
} }
...@@ -302,6 +303,7 @@ private: ...@@ -302,6 +303,7 @@ private:
std::vector<size_t>& suitableFaces = vertexToAdmissibleFaces[vertexID]; std::vector<size_t>& suitableFaces = vertexToAdmissibleFaces[vertexID];
if (suitableFaces.size()==0) { if (suitableFaces.size()==0) {
faultDofSet.erase(faultDofs.back());
faultDofs.pop_back(); faultDofs.pop_back();
faultFaces.pop_back(); faultFaces.pop_back();
vertexToAdmissibleFaces.erase(vertexID); vertexToAdmissibleFaces.erase(vertexID);
...@@ -327,6 +329,7 @@ private: ...@@ -327,6 +329,7 @@ private:
nextFaceDofs.erase(vertexID); nextFaceDofs.erase(vertexID);
faultDofs.push_back(*nextFaceDofs.begin()); faultDofs.push_back(*nextFaceDofs.begin());
faultDofSet.insert(*nextFaceDofs.begin());
faultFaces.push_back(nextFaceID); faultFaces.push_back(nextFaceID);
size_t nextVertexID = *(nextFaceDofs.begin()); size_t nextVertexID = *(nextFaceDofs.begin());
...@@ -352,19 +355,19 @@ private: ...@@ -352,19 +355,19 @@ private:
if (dir>0) { if (dir>0) {
while (lastDof<0) { while (lastDof<0) {
candidatIDs[dim]++;
auto dof = IDsToDof_[candidatIDs]; auto dof = IDsToDof_[candidatIDs];
if (separatingDofs.count(dof)) { if (separatingDofs.count(dof) or boundaryDofs_.count(dof)) {
lastDof = dof; lastDof = dof;
} }
candidatIDs[dim]++;
} }
} else { } else {
while (lastDof<0) { while (lastDof<0) {
candidatIDs[dim]--;
auto dof = IDsToDof_[candidatIDs]; auto dof = IDsToDof_[candidatIDs];
if (separatingDofs.count(dof)) { if (separatingDofs.count(dof) or boundaryDofs_.count(dof)) {
lastDof = dof; lastDof = dof;
} }
candidatIDs[dim]--;
} }
} }
...@@ -404,16 +407,16 @@ private: ...@@ -404,16 +407,16 @@ private:
void prolong(const MyRock& rock, MyRock& newRock) { void prolong(const MyRock& rock, MyRock& newRock) {
newRock.level = rock.level; 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(newLeft, newRight, newTop, newBottom, rock.center()); newRock.set(newLeft, newRight, newTop, newBottom, rock.center);
} }
bool randomSplit(const MyRock& rock) { bool randomSplit(const MyRock& rock) {
const auto& center = rock.center(); const auto& center = rock.center;
double res = 0.0; double res = 0.0;
if (center[0] < center[1]) { if (center[0] < center[1]) {
...@@ -433,16 +436,16 @@ private: ...@@ -433,16 +436,16 @@ private:
Rock newRock; Rock newRock;
prolong(rock, newRock); prolong(rock, newRock);
bool toBeSplit = (rock.level() == coarseLevelFactory_.level()) 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 ID 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] - vertexIDs_[newRock.left][0], size_t xCorridor = 1.0/4 * std::min(centerIDs[0] - vertexIDs_[newRock.left][0],
vertexIDs_[newRock.right][0] - centerIDs[0]) + 1; vertexIDs_[newRock.right][0] - centerIDs[0]);
size_t yCorridor = 1.0/2 * std::min(centerIDs[1] - vertexIDs_[newRock.bottom][1], size_t yCorridor = 1.0/4 * std::min(centerIDs[1] - vertexIDs_[newRock.bottom][1],
vertexIDs_[newRock.top][1] - centerIDs[1]) + 1; vertexIDs_[newRock.top][1] - centerIDs[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_);
...@@ -468,30 +471,30 @@ private: ...@@ -468,30 +471,30 @@ private:
MyRock rock00; MyRock rock00;
auto center00 = 1.0/2*(right + bottom); auto center00 = 1.0/2*(right + bottom);
createRock(rock00, center00, separatingDofs, x2Fault.getInterfaceDofs(), 0, y1Fault.getInterfaceDofs(), 0); createRock(rock00, center00, separatingDofs, x2Fault->getInterfaceDofs(), 0, y1Fault->getInterfaceDofs(), 0);
rocks_.push_back(rock00); rocks_.push_back(rock00);
MyRock rock01; MyRock rock01;
auto center01 = 1.0/2*(left + bottom); auto center01 = 1.0/2*(left + bottom);
createRock(rock01, center01, separatingDofs, x1Fault.getInterfaceDofs(), 0, y1Fault.getInterfaceDofs(), 1); createRock(rock01, center01, separatingDofs, x1Fault->getInterfaceDofs(), 0, y1Fault->getInterfaceDofs(), 1);
rocks_.push_back(rock01); rocks_.push_back(rock01);
MyRock rock10; MyRock rock10;
auto center10 = 1.0/2*(right + top); auto center10 = 1.0/2*(right + top);
createRock(rock10, center10, separatingDofs, x2Fault.getInterfaceDofs(), 1, y2Fault.getInterfaceDofs(), 0); createRock(rock10, center10, separatingDofs, x2Fault->getInterfaceDofs(), 1, y2Fault->getInterfaceDofs(), 0);
rocks_.push_back(rock10); rocks_.push_back(rock10);
MyRock rock11; MyRock rock11;
auto center11 = 1.0/2*(left + top); auto center11 = 1.0/2*(left + top);
createRock(rock11, center11, separatingDofs, x1Fault.getInterfaceDofs(), 1, y2Fault.getInterfaceDofs(), 1); createRock(rock11, center11, separatingDofs, x1Fault->getInterfaceDofs(), 1, y2Fault->getInterfaceDofs(), 1);
rocks_.push_back(rock11); rocks_.push_back(rock11);
} }
} }
public: public:
LevelRockFaultFactory(const std::shared_ptr<GridType> grid, const int level, const ctype resolution, LevelRockFaultFactory(const std::shared_ptr<GridType> grid, const int level, const ctype resolution,
const LevelRockFaultFactory& coarseLevelFactory, const LevelRockFaultFactory* coarseLevelFactory,
const double splittingThreshold = 0.0, const double maxAngle = 2) : const double splittingThreshold = 1.0, const double maxAngle = 2) :
grid_(grid), grid_(grid),
level_(level), level_(level),
resolution_(resolution), resolution_(resolution),
...@@ -500,14 +503,13 @@ private: ...@@ -500,14 +503,13 @@ private:
maxAngle_(maxAngle), maxAngle_(maxAngle),
indexSet_(gridView_.indexSet()), indexSet_(gridView_.indexSet()),
faceMapper_(gridView_), faceMapper_(gridView_),
coarseLevelFactory_(coarseLevelFactory) coarseLevelFactory_(*coarseLevelFactory)
{ {
// init faces_, vertexPositions_, vertexToFaces_ // init faces_, vertexPositions_, vertexToFaces_
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)); vertexIDs_.resize(gridView_.size(dim));
coarseToLevelVertex_.resize(coarseLevelFactory_.gridView().size(dim));
std::vector<bool> faceHandled(faceMapper_.size(), false); std::vector<bool> faceHandled(faceMapper_.size(), false);
...@@ -517,8 +519,11 @@ private: ...@@ -517,8 +519,11 @@ private:
const int faceID = faceMapper_.subIndex(elem, isect.indexInInside(), 1); const int faceID = faceMapper_.subIndex(elem, isect.indexInInside(), 1);
if (isect.boundary()) if (isect.boundary()) {
continue; std::set<size_t> intersectionDofs;
computeIntersectionDofs(isect, intersectionDofs);
boundaryDofs_.insert(intersectionDofs.begin(), intersectionDofs.end());
}
if (faceHandled[faceID]) if (faceHandled[faceID])
continue; continue;
...@@ -545,12 +550,15 @@ private: ...@@ -545,12 +550,15 @@ private:
} }
} }
for (const auto& vertex : vertices(gridView_)) { if (coarseLevelFactory!=nullptr) {
size_t coarseVertexID; const auto& coarseVertexCoords = coarseLevelFactory_.vertexCoords();
bool isCoarseVertex = coarseLevelFactory_.indexSet().index(vertex, coarseVertexID); coarseToLevelVertex_.resize(coarseVertexCoords.size());
for (size_t i=0; i<coarseVertexCoords.size(); i++) {
const auto& vertex = coarseVertexCoords[i];
ID id = {computeID(vertex, 0), computeID(vertex, 1)};
if (isCoarseVertex) { coarseToLevelVertex_[i] = IDsToDof_[id];
coarseToLevelVertex_[coarseVertexID] = indexSet_.index(vertex);
} }
} }
} }
...@@ -598,12 +606,17 @@ private: ...@@ -598,12 +606,17 @@ private:
const auto& faults() const { const auto& faults() const {
return faults_; return faults_;
} }
const auto& vertexCoords() const {
return vertexPositions_;
}
}; };
template <class GridType> template <class GridType>
class InitLevelRockFaultFactory : public LevelRockFaultFactory<GridType> { class InitLevelRockFaultFactory : public LevelRockFaultFactory<GridType> {
private: protected:
using Base = LevelRockFaultFactory<GridType>; using Base = LevelRockFaultFactory<GridType>;
public: public:
InitLevelRockFaultFactory(const std::shared_ptr<GridType> grid, const int level, const typename Base::ctype resolution) : InitLevelRockFaultFactory(const std::shared_ptr<GridType> grid, const int level, const typename Base::ctype resolution) :
Base(grid, level, resolution, nullptr) {} Base(grid, level, resolution, nullptr) {}
...@@ -634,6 +647,8 @@ template <class GridType> ...@@ -634,6 +647,8 @@ template <class GridType>
class RockFaultFactory { class RockFaultFactory {
private: private:
using GV = typename GridType::LevelGridView;
const int coarseResolution_; const int coarseResolution_;
const size_t maxLevel_; const size_t maxLevel_;
const int coarseGridN_; const int coarseGridN_;
...@@ -696,22 +711,23 @@ public: ...@@ -696,22 +711,23 @@ public:
// init level 0 rockFaultFactory // init level 0 rockFaultFactory
levelResolutions_[0] = std::pow(2, coarseResolution_); levelResolutions_[0] = std::pow(2, coarseResolution_);
InitLevelRockFaultFactory<GridType> initFactory(grid_, 0, levelResolutions_[0]);
std::set<size_t> boundaryDofs; std::set<size_t> boundaryDofs;
BoundaryIterator<GV> bIt(levelRockFaultFactories_[0].gridView(), BoundaryIterator<GV>::begin); BoundaryIterator<GV> bIt(initFactory.gridView(), BoundaryIterator<GV>::begin);
BoundaryIterator<GV> bEnd(levelRockFaultFactories_[0].gridView(), BoundaryIterator<GV>::end); BoundaryIterator<GV> bEnd(initFactory.gridView(), BoundaryIterator<GV>::end);
for(; bIt!=bEnd; ++bIt) { for(; bIt!=bEnd; ++bIt) {
std::set<size_t> intersectionDofs; std::set<size_t> intersectionDofs;
levelRockFaultFactories_[0].computeIntersectionDofs(*bIt, intersectionDofs); initFactory.computeIntersectionDofs(*bIt, intersectionDofs);
boundaryDofs.insert(intersectionDofs.begin(), intersectionDofs.end()); boundaryDofs.insert(intersectionDofs.begin(), intersectionDofs.end());
} }
InitLevelRockFaultFactory initFactory(grid_, 0, levelResolutions_[0]);
initFactory.build(boundaryDofs); initFactory.build(boundaryDofs);
levelRockFaultFactories_[0] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, 0, levelResolutions_[0], initFactory, 1.0); levelRockFaultFactories_[0] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, 0, levelResolutions_[0], &initFactory, 0.0);
levelRockFaultFactories_[0]->build(boundaryDofs); levelRockFaultFactories_[0]->build(boundaryDofs);
const auto& faults = levelRockFaultFactories_[0]->faults(); auto faults = levelRockFaultFactories_[0]->faults();
for (size_t j=0; j<faults.size(); j++) { for (size_t j=0; j<faults.size(); j++) {
interfaceNetwork_->addInterface(faults[j]); interfaceNetwork_->addInterface(faults[j]);
} }
...@@ -721,7 +737,7 @@ public: ...@@ -721,7 +737,7 @@ public:
levelResolutions_[i] = std::pow(2, coarseResolution_+i); levelResolutions_[i] = std::pow(2, coarseResolution_+i);
//generate faults on level //generate faults on level
levelRockFaultFactories_[i] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, i, levelResolutions_[i], *levelRockFaultFactories_[i-1], (i==1)*0.5); levelRockFaultFactories_[i] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, i, levelResolutions_[i], levelRockFaultFactories_[i-1].get(), (i==1)*0.5 + (i>1)*1.0);
levelRockFaultFactories_[i]->build(interfaceNetwork_->getInterfaceNetworkDofs(i)); levelRockFaultFactories_[i]->build(interfaceNetwork_->getInterfaceNetworkDofs(i));
faults = levelRockFaultFactories_[i]->faults(); faults = levelRockFaultFactories_[i]->faults();
......
src/cantorfaultnetworks/results/sparse/plot.tex
+ 18
− 21
View file @ b0c4be5c
...@@ -4,25 +4,22 @@ ...@@ -4,25 +4,22 @@
\begin{document} \begin{document}
\def\scale{6}
\begin{figure} \begin{minipage}[t]{0.45\linewidth}
\def\scale{5} \flushleft
\input{levelinterfacenetwork_0.tikz} \input{levelinterfacenetwork_0.tikz}\\
\end{figure} \vspace{2em}
\input{levelinterfacenetwork_1.tikz}\\
\begin{figure} \vspace{2em}
\def\scale{5} \input{levelinterfacenetwork_2.tikz}\\
\input{levelinterfacenetwork_2.tikz} \vfill
\end{figure} \end{minipage}
\hspace{3em}
\begin{figure} \begin{minipage}[t]{0.45\linewidth}
\def\scale{5} \flushright
\input{levelinterfacenetwork_4.tikz} \input{levelinterfacenetwork_3.tikz}\\
\end{figure} \vspace{2em}
%\input{levelinterfacenetwork_4.tikz}
% \begin{figure} \vfill
% \def\scale{5} \end{minipage}
% \input{levelinterfacenetwork_6.tikz}
% \end{figure}
\end{document} \end{document}
\ No newline at end of file
src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
+ 6
− 2
View file @ b0c4be5c
...@@ -152,11 +152,15 @@ int main(int argc, char** argv) { try ...@@ -152,11 +152,15 @@ int main(int argc, char** argv) { try
const int exactLevelIdx = 2*(maxCantorLevel - minCantorLevel); const int exactLevelIdx = 2*(maxCantorLevel - minCantorLevel);
// build multilevel cantor fault network // build multilevel cantor fault network
SparseCantorFaultFactory<GridType> faultFactory(minCantorLevel, maxCantorLevel); //SparseCantorFaultFactory<GridType> faultFactory(minCantorLevel, maxCantorLevel);
//const auto& interfaceNetwork = faultFactory.interfaceNetwork();
std::srand(5);
RockFaultFactory<GridType> faultFactory(5, 1);
const auto& interfaceNetwork = faultFactory.interfaceNetwork(); const auto& interfaceNetwork = faultFactory.interfaceNetwork();
if (problemCount==0) { if (problemCount==0) {
std::vector<int> writeLevels {0, 2}; std::vector<int> writeLevels {0, 1, 2};
for (size_t i=0; i<writeLevels.size(); i++) { for (size_t i=0; i<writeLevels.size(); i++) {
int writeLevel = writeLevels[i]; int writeLevel = writeLevels[i];
bool writeGrid = !(writeLevel==8); bool writeGrid = !(writeLevel==8);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment