From 10709c3816b6d0ef444f3092043a97786f73de27 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Fri, 12 Jun 2020 14:25:12 +0200
Subject: [PATCH] .
---
.../faultfactories/rockfaultfactory.hh | 79 +++++++++++++------
1 file changed, 56 insertions(+), 23 deletions(-)
diff --git a/dune/faultnetworks/faultfactories/rockfaultfactory.hh b/dune/faultnetworks/faultfactories/rockfaultfactory.hh
index e48d563..81404e8 100644
--- a/dune/faultnetworks/faultfactories/rockfaultfactory.hh
+++ b/dune/faultnetworks/faultfactories/rockfaultfactory.hh
@@ -121,7 +121,7 @@ private:
// check if intersection has separating dofs or other fault dofs
std::set<size_t> intersectionDofs;
- computeIntersectionDofs(indexSet_, intersection, intersectionDofs);
+ computeIntersectionDofs(intersection, intersectionDofs);
intersectionDofs.erase(vertexID);
const size_t otherDim = (dim + 1) % 2;
@@ -151,21 +151,6 @@ private:
return (int) (vertex[dim]*resolution_);
}
- void computeIntersectionDofs(const Intersection& intersection, std::set<size_t>& intersectionDofs) const {
- intersectionDofs.clear();
-
- // loop over all vertices of the intersection
- const auto& insideElement = intersection.inside();
- const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
-
- for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
- size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
- size_t globalIdx = indexSet_.subIndex(insideElement, idxInElement, dimElement);
-
- intersectionDofs.insert(globalIdx);
- }
- }
-
bool isTargetBoundaryVertex(size_t vertexIdx, const std::array<size_t, 2>& targetIDs) const {
const auto& vertex = vertexPositions_[vertexIdx];
const std::array<size_t, 2> vertexIDs = {computeID(vertex, 0), computeID(vertex, 1)};
@@ -338,7 +323,7 @@ private:
suitableFaces.erase(suitableFaces.begin()+randSelector);
std::set<size_t> nextFaceDofs;
- computeIntersectionDofs(indexSet_, faces_[nextFaceID], nextFaceDofs);
+ computeIntersectionDofs(faces_[nextFaceID], nextFaceDofs);
nextFaceDofs.erase(vertexID);
faultDofs.push_back(*nextFaceDofs.begin());
@@ -504,13 +489,18 @@ private:
}
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 double splittingThreshold = 0.0, const double maxAngle = 2) :
grid_(grid),
level_(level),
resolution_(resolution),
gridView_(grid_->levelGridView(level_)),
+ splittingThreshold_(splittingThreshold),
+ maxAngle_(maxAngle),
indexSet_(gridView_.indexSet()),
- faceMapper_(gridView_)
+ faceMapper_(gridView_),
+ coarseLevelFactory_(coarseLevelFactory)
{
// init faces_, vertexPositions_, vertexToFaces_
faces_.resize(faceMapper_.size());
@@ -574,6 +564,25 @@ private:
}
}
+ void computeIntersectionDofs(const Intersection& intersection, std::set<size_t>& intersectionDofs) const {
+ intersectionDofs.clear();
+
+ // loop over all vertices of the intersection
+ const auto& insideElement = intersection.inside();
+ const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
+
+ for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
+ size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
+ size_t globalIdx = indexSet_.subIndex(insideElement, idxInElement, dimElement);
+
+ intersectionDofs.insert(globalIdx);
+ }
+ }
+
+ const auto& gridView() const {
+ return gridView_;
+ }
+
const auto& indexSet() const {
return indexSet_;
}
@@ -585,6 +594,10 @@ private:
auto level() const {
return level_;
}
+
+ const auto& faults() const {
+ return faults_;
+ }
};
template <class GridType>
@@ -650,15 +663,35 @@ public:
// init interface network
interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);
- for (size_t i=0; i<=maxLevel_; i++) {
- std::vector<std::shared_ptr<FaultInterface<GV>>> faults;
+ // init level 0 rockFaultFactory
+ levelResolutions_[0] = std::pow(2, coarseResolution_);
+ levelRockFaultFactories_[0] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, 0, levelResolutions_[0], levelRockFaultFactories_[i-1], 1.0);
+
+ std::set<size_t> boundaryDofs;
+ BoundaryIterator<GV> bIt(levelRockFaultFactories_[0].gridView(), BoundaryIterator<GV>::begin);
+ BoundaryIterator<GV> bEnd(levelRockFaultFactories_[0].gridView(), BoundaryIterator<GV>::end);
+ for(; bIt!=bEnd; ++bIt) {
+ std::set<size_t> intersectionDofs;
+ levelRockFaultFactories_[0].computeIntersectionDofs(*bIt, intersectionDofs);
+ boundaryDofs.insert(intersectionDofs.begin(), intersectionDofs.end());
+ }
+
+ levelRockFaultFactories_[0]->build(boundaryDofs);
+
+ const auto& faults = levelRockFaultFactories_[0]->faults();
+ for (size_t j=0; j<faults.size(); j++) {
+ interfaceNetwork_->addInterface(faults[j]);
+ }
+ interfaceNetwork_->prolongLevel(0, 1);
+ for (size_t i=1; i<=maxLevel_; i++) {
levelResolutions_[i] = std::pow(2, coarseResolution_+i);
//generate faults on level
- levelRockFaultFactories_[i] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, i, levelResolutions_[i]);
- levelRockFaultFactories_[i]->build(faults, maxAngle);
+ levelRockFaultFactories_[i] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, i, levelResolutions_[i], levelRockFaultFactories_[i-1], (i==1)*0.5);
+ levelRockFaultFactories_[i]->build(interfaceNetwork_->getInterfaceNetworkDofs(i));
+ faults = levelRockFaultFactories_[i]->faults();
for (size_t j=0; j<faults.size(); j++) {
interfaceNetwork_->addInterface(faults[j]);
}
--
GitLab