.

dune/faultnetworks/faultfactories/rockfaultfactory.hh

@@ -378,7 +378,7 @@ private:
         rock.level = level_;
         rock.center = center;
 
-        const ID centerIDs = {computeID(center(), 0), computeID(center(), 1)};
+        const ID centerIDs = {computeID(center, 0), computeID(center, 1)};
 
         size_t left, right, top, bottom = 0;
 
@@ -600,6 +600,36 @@ private:
         }
 };
 
+template <class GridType>
+class InitLevelRockFaultFactory : public LevelRockFaultFactory<GridType> {
+private:
+    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) {}
+
+    void build(const std::set<size_t>& boundaryDofs) {
+        typename Base::MyRock rock;
+
+        rock.level = this->level_;
+        rock.center[0] = 0.5;
+        rock.center[1] = 0.5;
+
+        const typename Base::ID centerIDs = {this->computeID(rock.center, 0), this->computeID(rock.center, 1)};
+
+        size_t top = this->searchDof(centerIDs, boundaryDofs, 1, 1);
+        size_t bottom = this->searchDof(centerIDs, boundaryDofs, 1, -1);
+
+        size_t left = this->searchDof(centerIDs, boundaryDofs, 0, -1);
+        size_t right = this->searchDof(centerIDs, boundaryDofs, 0, 1);
+
+        rock.set(left, right, top, bottom, rock.center);
+
+        this->rocks_.push_back(rock);
+    }
+
+};
+
 template <class GridType>
 class RockFaultFactory {
 
@@ -665,7 +695,6 @@ public:
 
             // 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);
@@ -676,6 +705,10 @@ public:
                 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]->build(boundaryDofs);
 
             const auto& faults = levelRockFaultFactories_[0]->faults();
@@ -688,7 +721,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], (i==1)*0.5);
                 levelRockFaultFactories_[i]->build(interfaceNetwork_->getInterfaceNetworkDofs(i));
 
                 faults = levelRockFaultFactories_[i]->faults();