.

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]);
                 }