From b0c4be5c9fbd6572c5c048173d1cb5c68740d0ef Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Sun, 14 Jun 2020 17:36:46 +0200
Subject: [PATCH] .
---
.../faultfactories/rockfaultfactory.hh | 118 ++++++++++--------
.../results/sparse/plot.tex | 39 +++---
.../sparsecantorfaultnetwork.cc | 8 +-
3 files changed, 91 insertions(+), 74 deletions(-)
diff --git a/dune/faultnetworks/faultfactories/rockfaultfactory.hh b/dune/faultnetworks/faultfactories/rockfaultfactory.hh
index 7f5bdef..31e0f13 100644
--- a/dune/faultnetworks/faultfactories/rockfaultfactory.hh
+++ b/dune/faultnetworks/faultfactories/rockfaultfactory.hh
@@ -91,22 +91,20 @@ protected:
std::vector<ID> vertexIDs_;
std::map<ID, size_t> IDsToDof_;
+ std::set<size_t> boundaryDofs_;
+
std::vector<std::vector<size_t>> vertexToFaces_;
std::vector<int> coarseToLevelVertex_;
const LevelRockFaultFactory& coarseLevelFactory_;
std::vector<MyRock> rocks_;
- std::vector<std::shared_ptr<FaultInterface<GV>>>& faults_;
-
-private:
- template <typename T> int sgn(T val) {
- return (T(0) < val) - (val < T(0));
- }
+ std::vector<std::shared_ptr<FaultInterface<GV>>> faults_;
+protected:
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 std::set<size_t>& separatingIDs, const std::set<size_t>& faultDofs,
const Coords& direction, size_t dim) const {
const auto& intersection = faces_[faceID];
@@ -127,9 +125,10 @@ private:
const size_t otherDim = (dim + 1) % 2;
for (const auto& isDof : intersectionDofs) {
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)) {
return false;
}
@@ -250,17 +249,19 @@ private:
if (std::abs(direction[0]) < std::abs(direction[1])) {
dim = 1;
}
-
+ auto otherDim = (dim+1) % 2;
const ID centerIDs = {computeID(center, 0), computeID(center, 1)};
std::set<size_t> separatingIDs;
- auto faultSeedID = vertexIDs_[faultSeedIdx][dim];
+ auto faultSeedID = vertexIDs_[faultSeedIdx][otherDim];
separatingIDs.insert(faultSeedID + corridor);
separatingIDs.insert(faultSeedID - corridor);
std::vector<size_t> faultDofs(0);
+ std::set<size_t> faultDofSet;
faultDofs.push_back(faultSeedIdx);
+ faultDofSet.insert(faultSeedIdx);
std::vector<size_t> faultFaces;
@@ -283,7 +284,7 @@ private:
std::vector<size_t> admissibleFaces;
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]);
}
}
@@ -302,6 +303,7 @@ private:
std::vector<size_t>& suitableFaces = vertexToAdmissibleFaces[vertexID];
if (suitableFaces.size()==0) {
+ faultDofSet.erase(faultDofs.back());
faultDofs.pop_back();
faultFaces.pop_back();
vertexToAdmissibleFaces.erase(vertexID);
@@ -327,6 +329,7 @@ private:
nextFaceDofs.erase(vertexID);
faultDofs.push_back(*nextFaceDofs.begin());
+ faultDofSet.insert(*nextFaceDofs.begin());
faultFaces.push_back(nextFaceID);
size_t nextVertexID = *(nextFaceDofs.begin());
@@ -352,19 +355,19 @@ private:
if (dir>0) {
while (lastDof<0) {
- candidatIDs[dim]++;
auto dof = IDsToDof_[candidatIDs];
- if (separatingDofs.count(dof)) {
+ if (separatingDofs.count(dof) or boundaryDofs_.count(dof)) {
lastDof = dof;
}
+ candidatIDs[dim]++;
}
} else {
while (lastDof<0) {
- candidatIDs[dim]--;
auto dof = IDsToDof_[candidatIDs];
- if (separatingDofs.count(dof)) {
+ if (separatingDofs.count(dof) or boundaryDofs_.count(dof)) {
lastDof = dof;
}
+ candidatIDs[dim]--;
}
}
@@ -404,16 +407,16 @@ private:
void prolong(const MyRock& rock, MyRock& 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()];
+ auto newLeft = coarseToLevelVertex_[rock.left];
+ auto newRight = coarseToLevelVertex_[rock.right];
+ auto newTop = coarseToLevelVertex_[rock.top];
+ 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) {
- const auto& center = rock.center();
+ const auto& center = rock.center;
double res = 0.0;
if (center[0] < center[1]) {
@@ -433,16 +436,16 @@ private:
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) {
rocks_.push_back(newRock);
} else {
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],
- 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;
+ size_t xCorridor = 1.0/4 * std::min(centerIDs[0] - vertexIDs_[newRock.left][0],
+ vertexIDs_[newRock.right][0] - centerIDs[0]);
+ size_t yCorridor = 1.0/4 * std::min(centerIDs[1] - vertexIDs_[newRock.bottom][1],
+ vertexIDs_[newRock.top][1] - centerIDs[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_);
@@ -468,30 +471,30 @@ private:
MyRock rock00;
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);
MyRock rock01;
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);
MyRock rock10;
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);
MyRock rock11;
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);
}
}
public:
LevelRockFaultFactory(const std::shared_ptr<GridType> grid, const int level, const ctype resolution,
- const LevelRockFaultFactory& coarseLevelFactory,
- const double splittingThreshold = 0.0, const double maxAngle = 2) :
+ const LevelRockFaultFactory* coarseLevelFactory,
+ const double splittingThreshold = 1.0, const double maxAngle = 2) :
grid_(grid),
level_(level),
resolution_(resolution),
@@ -500,14 +503,13 @@ private:
maxAngle_(maxAngle),
indexSet_(gridView_.indexSet()),
faceMapper_(gridView_),
- coarseLevelFactory_(coarseLevelFactory)
+ coarseLevelFactory_(*coarseLevelFactory)
{
// init faces_, vertexPositions_, vertexToFaces_
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);
@@ -517,8 +519,11 @@ private:
const int faceID = faceMapper_.subIndex(elem, isect.indexInInside(), 1);
- if (isect.boundary())
- continue;
+ if (isect.boundary()) {
+ std::set<size_t> intersectionDofs;
+ computeIntersectionDofs(isect, intersectionDofs);
+ boundaryDofs_.insert(intersectionDofs.begin(), intersectionDofs.end());
+ }
if (faceHandled[faceID])
continue;
@@ -545,15 +550,18 @@ private:
}
}
- for (const auto& vertex : vertices(gridView_)) {
- size_t coarseVertexID;
- bool isCoarseVertex = coarseLevelFactory_.indexSet().index(vertex, coarseVertexID);
+ if (coarseLevelFactory!=nullptr) {
+ const auto& coarseVertexCoords = coarseLevelFactory_.vertexCoords();
+ coarseToLevelVertex_.resize(coarseVertexCoords.size());
- if (isCoarseVertex) {
- coarseToLevelVertex_[coarseVertexID] = indexSet_.index(vertex);
+ for (size_t i=0; i<coarseVertexCoords.size(); i++) {
+ const auto& vertex = coarseVertexCoords[i];
+ ID id = {computeID(vertex, 0), computeID(vertex, 1)};
+
+ coarseToLevelVertex_[i] = IDsToDof_[id];
}
}
- }
+ }
void build(const std::set<size_t>& separatingDofs) {
faults_.resize(0);
@@ -598,12 +606,17 @@ private:
const auto& faults() const {
return faults_;
}
+
+ const auto& vertexCoords() const {
+ return vertexPositions_;
+ }
};
template <class GridType>
class InitLevelRockFaultFactory : public LevelRockFaultFactory<GridType> {
-private:
+protected:
using Base = LevelRockFaultFactory<GridType>;
+
public:
InitLevelRockFaultFactory(const std::shared_ptr<GridType> grid, const int level, const typename Base::ctype resolution) :
Base(grid, level, resolution, nullptr) {}
@@ -634,6 +647,8 @@ template <class GridType>
class RockFaultFactory {
private:
+ using GV = typename GridType::LevelGridView;
+
const int coarseResolution_;
const size_t maxLevel_;
const int coarseGridN_;
@@ -696,22 +711,23 @@ public:
// init level 0 rockFaultFactory
levelResolutions_[0] = std::pow(2, coarseResolution_);
+ InitLevelRockFaultFactory<GridType> initFactory(grid_, 0, levelResolutions_[0]);
+
std::set<size_t> boundaryDofs;
- BoundaryIterator<GV> bIt(levelRockFaultFactories_[0].gridView(), BoundaryIterator<GV>::begin);
- BoundaryIterator<GV> bEnd(levelRockFaultFactories_[0].gridView(), BoundaryIterator<GV>::end);
+ BoundaryIterator<GV> bIt(initFactory.gridView(), BoundaryIterator<GV>::begin);
+ BoundaryIterator<GV> bEnd(initFactory.gridView(), BoundaryIterator<GV>::end);
for(; bIt!=bEnd; ++bIt) {
std::set<size_t> intersectionDofs;
- levelRockFaultFactories_[0].computeIntersectionDofs(*bIt, intersectionDofs);
+ initFactory.computeIntersectionDofs(*bIt, intersectionDofs);
boundaryDofs.insert(intersectionDofs.begin(), intersectionDofs.end());
}
- InitLevelRockFaultFactory initFactory(grid_, 0, levelResolutions_[0]);
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);
- const auto& faults = levelRockFaultFactories_[0]->faults();
+ auto faults = levelRockFaultFactories_[0]->faults();
for (size_t j=0; j<faults.size(); j++) {
interfaceNetwork_->addInterface(faults[j]);
}
@@ -721,7 +737,7 @@ public:
levelResolutions_[i] = std::pow(2, coarseResolution_+i);
//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));
faults = levelRockFaultFactories_[i]->faults();
diff --git a/src/cantorfaultnetworks/results/sparse/plot.tex b/src/cantorfaultnetworks/results/sparse/plot.tex
index 5c038a7..6633d7c 100644
--- a/src/cantorfaultnetworks/results/sparse/plot.tex
+++ b/src/cantorfaultnetworks/results/sparse/plot.tex
@@ -4,25 +4,22 @@
\begin{document}
-
- \begin{figure}
- \def\scale{5}
- \input{levelinterfacenetwork_0.tikz}
- \end{figure}
-
- \begin{figure}
- \def\scale{5}
- \input{levelinterfacenetwork_2.tikz}
- \end{figure}
-
- \begin{figure}
- \def\scale{5}
- \input{levelinterfacenetwork_4.tikz}
- \end{figure}
-
-% \begin{figure}
-% \def\scale{5}
-% \input{levelinterfacenetwork_6.tikz}
-% \end{figure}
-
+ \def\scale{6}
+ \begin{minipage}[t]{0.45\linewidth}
+ \flushleft
+ \input{levelinterfacenetwork_0.tikz}\\
+ \vspace{2em}
+ \input{levelinterfacenetwork_1.tikz}\\
+ \vspace{2em}
+ \input{levelinterfacenetwork_2.tikz}\\
+ \vfill
+ \end{minipage}
+ \hspace{3em}
+ \begin{minipage}[t]{0.45\linewidth}
+ \flushright
+ \input{levelinterfacenetwork_3.tikz}\\
+ \vspace{2em}
+ %\input{levelinterfacenetwork_4.tikz}
+ \vfill
+ \end{minipage}
\end{document}
\ No newline at end of file
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
index e32efc7..9297708 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
@@ -152,11 +152,15 @@ int main(int argc, char** argv) { try
const int exactLevelIdx = 2*(maxCantorLevel - minCantorLevel);
// 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();
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++) {
int writeLevel = writeLevels[i];
bool writeGrid = !(writeLevel==8);
--
GitLab