.

dune/faultnetworks/faultfactories/rockfaultfactory.hh

@@ -127,8 +127,11 @@ protected:
                 const auto& vertex = vertexPositions_[isDof];
                 ID IDs = {computeID(vertex, 0), computeID(vertex, 1)};
 
+                if ((centerIDs[0]==IDs[0]) and (centerIDs[1]==IDs[1]))
+                        return true;
+
                 auto val = direction[dim];
-                bool centerPassed = ((0 < val) - (val < 0)>0) ? (centerIDs[dim]<IDs[dim]) : (centerIDs[dim]>IDs[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;
                 }
@@ -147,6 +150,7 @@ protected:
 
         //works only in 2D, vertex(1) = vertex(0) - yid, intersection of line (given by vertex and derivative 1) with y axis
         int computeID(const Coords& vertex, int dim) const {
+            auto v = vertex;
             return (int) (vertex[dim]*resolution_);
         }
 
@@ -231,7 +235,7 @@ protected:
         }*/
 
 
-        void generateFault(std::shared_ptr<FaultInterface<GV>> fault, const size_t faultSeedIdx, const Coords& center,
+        bool generateFault(std::shared_ptr<FaultInterface<GV>> fault, const size_t faultSeedIdx, const Coords& center,
                            const size_t corridor,
                            const std::set<size_t>& separatingDofs) {
             //std::cout << "LevelGeoFaultFactory::generateFault() " << std::endl;
@@ -263,7 +267,7 @@ protected:
             faultDofs.push_back(faultSeedIdx);
             faultDofSet.insert(faultSeedIdx);
 
-            std::vector<size_t> faultFaces;
+            std::vector<size_t> faultFaces(0);
 
             std::map<size_t, std::vector<size_t>> vertexToAdmissibleFaces;
 
@@ -303,13 +307,18 @@ protected:
                 std::vector<size_t>& suitableFaces = vertexToAdmissibleFaces[vertexID];
 
                 if (suitableFaces.size()==0) {
-                    faultDofSet.erase(faultDofs.back());
-                    faultDofs.pop_back();
-                    faultFaces.pop_back();
-                    vertexToAdmissibleFaces.erase(vertexID);
+                    if (faultDofs.size()>0) {
+                        faultDofSet.erase(faultDofs.back());
+                        faultDofs.pop_back();
+                        faultFaces.pop_back();
+                        vertexToAdmissibleFaces.erase(vertexID);
+
+                        deadendDofs.insert(vertexID);
+                    }
 
-                    deadendDofs.insert(vertexID);
-                    vertexQueue.push(faultDofs.back());
+                    if (faultDofs.size()>0) {
+                        vertexQueue.push(faultDofs.back());
+                    }
                 } else {
 
                     // generate random number from (0,1)
@@ -337,16 +346,13 @@ protected:
                 }
             }
 
-            if (!success) {
-                std::cout << "Generating a fault failed: Unable to reach target boundary! This should not happend!" << std::endl;
-                DUNE_THROW(Dune::Exception, "Generating a fault failed: Unable to reach target boundary!");
-            }
-
-            for (size_t i=0; i<faultFaces.size(); i++) {
-                fault->addFace(faces_[faultFaces[i]]);
+            if (success) {
+                for (size_t i=0; i<faultFaces.size(); i++) {
+                    fault->addFace(faces_[faultFaces[i]]);
+                }
             }
 
-            //std::cout << "------------------------------------- " << std::endl << std::endl;
+            return success;
         }
 
     auto searchDof(const ID& IDs, const std::set<size_t>& separatingDofs, size_t dim, int dir) {
@@ -374,12 +380,12 @@ protected:
         return lastDof;
     }
 
+
     void createRock(MyRock& rock, const Coords& center, const std::set<size_t>& separatingDofs,
                     const std::set<size_t>& xFaultDofs, bool xFaultBottom,
                     const std::set<size_t>& yFaultDofs, bool yFaultRight) {
 
         rock.level = level_;
-        rock.center = center;
 
         const ID centerIDs = {computeID(center, 0), computeID(center, 1)};
 
@@ -432,69 +438,98 @@ protected:
         return  100.0*res > prob;
     }
 
+    void appendFaultDofs(const std::shared_ptr<FaultInterface<GV>>& fault, const size_t centerIdx, std::set<size_t>& dofs) {
+        const auto& faultDofs = fault->getInterfaceDofs();
+        dofs.insert(faultDofs.begin(), faultDofs.end());
+        dofs.erase(centerIdx);
+    }
+
     void split(const MyRock& rock, const std::set<size_t>& separatingDofs) {
+        auto forbiddenDofs = separatingDofs;
+
         Rock newRock;
         prolong(rock, newRock);
 
+        std::set<int> set;
+        set.insert(newRock.left);
+        set.insert(newRock.right);
+        set.insert(newRock.top);
+        set.insert(newRock.bottom);
+        bool canSplit = (set.size()==4);
+
         bool toBeSplit = (rock.level == coarseLevelFactory_.level()) and randomSplit(rock);
-        if (!toBeSplit) {
+        if (!toBeSplit or !canSplit) {
             rocks_.push_back(newRock);
         } else {
             const ID centerIDs = {computeID(newRock.center, 0), computeID(newRock.center, 1)};
+            const size_t centerIdx = IDsToDof_[centerIDs];
 
-            size_t xCorridor = 1.0/4 * std::min(centerIDs[0] - vertexIDs_[newRock.left][0],
+            size_t xCorridor = 1.0/2 * 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],
+            size_t yCorridor = 1.0/2 * std::min(centerIDs[1] - vertexIDs_[newRock.bottom][1],
                     vertexIDs_[newRock.top][1] - centerIDs[1]);
 
+            auto& l = vertexIDs_[newRock.left];
+            auto& r = vertexIDs_[newRock.right];
+            auto& t = vertexIDs_[newRock.top];
+            auto& b = vertexIDs_[newRock.bottom];
+
             // 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_);
-            generateFault(x1Fault, newRock.left, newRock.center, yCorridor, separatingDofs);
-            faults_.push_back(x1Fault);
+            auto x1 = generateFault(x1Fault, newRock.left, newRock.center, yCorridor, forbiddenDofs);
+            appendFaultDofs(x1Fault, centerIdx, forbiddenDofs);
 
             std::shared_ptr<FaultInterface<GV>> x2Fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
-            generateFault(x2Fault, newRock.right, newRock.center, yCorridor, separatingDofs);
-            faults_.push_back(x2Fault);
+            auto x2 = generateFault(x2Fault, newRock.right, newRock.center, yCorridor, forbiddenDofs);
+            appendFaultDofs(x2Fault, centerIdx, forbiddenDofs);
 
             std::shared_ptr<FaultInterface<GV>> y1Fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
-            generateFault(y1Fault, newRock.bottom, newRock.center, xCorridor, separatingDofs);
-            faults_.push_back(y1Fault);
+            auto y1 = generateFault(y1Fault, newRock.bottom, newRock.center, xCorridor, forbiddenDofs);
+            appendFaultDofs(y1Fault, centerIdx, forbiddenDofs);
 
             std::shared_ptr<FaultInterface<GV>> y2Fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
-            generateFault(y2Fault, newRock.top, newRock.center, xCorridor, separatingDofs);
-            faults_.push_back(y2Fault);
-
-            auto left = vertexPositions_[newRock.left];
-            auto right = vertexPositions_[newRock.right];
-            auto top = vertexPositions_[newRock.top];
-            auto bottom = vertexPositions_[newRock.bottom];
-
-            MyRock rock00;
-            auto center00 = 1.0/2*(right + bottom);
-            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);
-            rocks_.push_back(rock01);
-
-            MyRock rock10;
-            auto center10 = 1.0/2*(right + top);
-            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);
-            rocks_.push_back(rock11);
+            auto y2 = generateFault(y2Fault, newRock.top, newRock.center, xCorridor, forbiddenDofs);
+
+            bool generatingFaultsSuccessful = x1 and x2 and y1 and y2;
+
+            if (generatingFaultsSuccessful) {
+                faults_.push_back(x1Fault);
+                faults_.push_back(x2Fault);
+                faults_.push_back(y1Fault);
+                faults_.push_back(y2Fault);
+
+                auto left = vertexPositions_[newRock.left];
+                auto right = vertexPositions_[newRock.right];
+                auto top = vertexPositions_[newRock.top];
+                auto bottom = vertexPositions_[newRock.bottom];
+
+                MyRock rock00;
+                auto center00 = 1.0/2*(right + bottom);
+                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);
+                rocks_.push_back(rock01);
+
+                MyRock rock10;
+                auto center10 = 1.0/2*(right + top);
+                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);
+                rocks_.push_back(rock11);
+            }
         }
     }
 
     public:
         LevelRockFaultFactory(const std::shared_ptr<GridType> grid, const int level, const ctype resolution,
-                              const LevelRockFaultFactory* coarseLevelFactory,
-                              const double splittingThreshold = 1.0, const double maxAngle = 2) :
+                              const std::shared_ptr<LevelRockFaultFactory> coarseLevelFactory,
+                              const double splittingThreshold = 1.0, const double maxAngle = 2.0) :
             grid_(grid),
             level_(level),
             resolution_(resolution),
@@ -519,18 +554,19 @@ protected:
 
                         const int faceID = faceMapper_.subIndex(elem, isect.indexInInside(), 1);
 
-                        if (isect.boundary()) {
-                            std::set<size_t> intersectionDofs;
-                            computeIntersectionDofs(isect, intersectionDofs);
-                            boundaryDofs_.insert(intersectionDofs.begin(), intersectionDofs.end());
-                        }
-
                         if (faceHandled[faceID])
                             continue;
 
                         faceHandled[faceID] = true;
                         faces_[faceID] = isect;
 
+                        if (isect.boundary()) {
+                            std::set<size_t> intersectionDofs;
+                            computeIntersectionDofs(isect, intersectionDofs);
+                            boundaryDofs_.insert(intersectionDofs.begin(), intersectionDofs.end());
+                            continue;
+                        }
+
                         const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(elem.type());
 
                         for (int i=0; i<refElement.size(isect.indexInInside(), 1, dimElement); i++) {
@@ -711,20 +747,20 @@ public:
             // init level 0 rockFaultFactory
             levelResolutions_[0] = std::pow(2, coarseResolution_);
 
-            InitLevelRockFaultFactory<GridType> initFactory(grid_, 0, levelResolutions_[0]);
+            auto initFactory = std::make_shared<InitLevelRockFaultFactory<GridType>>(grid_, 0, levelResolutions_[0]);
 
             std::set<size_t> boundaryDofs;
-            BoundaryIterator<GV> bIt(initFactory.gridView(), BoundaryIterator<GV>::begin);
-            BoundaryIterator<GV> bEnd(initFactory.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;
-                initFactory.computeIntersectionDofs(*bIt, intersectionDofs);
+                initFactory->computeIntersectionDofs(*bIt, intersectionDofs);
                 boundaryDofs.insert(intersectionDofs.begin(), intersectionDofs.end());
             }
 
-            initFactory.build(boundaryDofs);
+            initFactory->build(boundaryDofs);
 
-            levelRockFaultFactories_[0] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, 0, levelResolutions_[0], &initFactory, 0.0);
+            levelRockFaultFactories_[0] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, 0, levelResolutions_[0], initFactory, 0.0, maxAngle);
             levelRockFaultFactories_[0]->build(boundaryDofs);
 
             auto faults = levelRockFaultFactories_[0]->faults();
@@ -737,7 +773,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].get(), (i==1)*0.5 + (i>1)*1.0);
+                levelRockFaultFactories_[i] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, i, levelResolutions_[i], levelRockFaultFactories_[i-1], (i==1)*0.5 + (i>1)*1.0, maxAngle);
                 levelRockFaultFactories_[i]->build(interfaceNetwork_->getInterfaceNetworkDofs(i));
 
                 faults = levelRockFaultFactories_[i]->faults();

dune/faultnetworks/faultp1nodalbasis.hh

@@ -17,6 +17,8 @@
 
 #include <dune/faultnetworks/levelinterfacenetwork.hh>
 
+#include <dune/faultnetworks/utils/debugutils.hh>
+
 template <class GV, class RT=double>
 class FaultP1NodalBasis :
     public P1NodalBasis<

src/cantorfaultnetworks/results/sparse/plot.tex

@@ -17,7 +17,7 @@
 	\hspace{3em}
 	\begin{minipage}[t]{0.45\linewidth}	
 			\flushright
-			\input{levelinterfacenetwork_3.tikz}\\
+			%\input{levelinterfacenetwork_3.tikz}\\
 			\vspace{2em}
 			%\input{levelinterfacenetwork_4.tikz}
 			\vfill

src/cantorfaultnetworks/sparsecantorfaultnetwork.cc

@@ -59,7 +59,6 @@
 #include <dune/faultnetworks/dgmgtransfer.hh>
 
 #include <dune/faultnetworks/faultfactories/sparsecantorfaultfactory.hh>
-#include <dune/faultnetworks/faultfactories/rockfaultfactory.hh>
 
 #include <dune/solvers/solvers/loopsolver.hh>
 
@@ -152,11 +151,7 @@ int main(int argc, char** argv) { try
     const int exactLevelIdx = 2*(maxCantorLevel - minCantorLevel);
 
     // build multilevel cantor fault network
-    //SparseCantorFaultFactory<GridType> faultFactory(minCantorLevel, maxCantorLevel);
-    //const auto& interfaceNetwork = faultFactory.interfaceNetwork();
-
-    std::srand(5);
-    RockFaultFactory<GridType> faultFactory(5, 1);
+    SparseCantorFaultFactory<GridType> faultFactory(minCantorLevel, maxCantorLevel);
     const auto& interfaceNetwork = faultFactory.interfaceNetwork();
 
     if (problemCount==0) {

src/geofaultnetworks/CMakeLists.txt

 add_executable("geofaultnetwork" geofaultnetwork.cc)
 target_link_dune_default_libraries("geofaultnetwork")
 
+add_executable("rockfaultnetwork" rockfaultnetwork.cc)
+target_link_dune_default_libraries("rockfaultnetwork")
+
 add_custom_target(geofaultnetworks_srcs SOURCES
 geofaultnetwork.parset
 )
 
-
+add_custom_target(rockfaultnetworks_srcs SOURCES
+rockfaultnetwork.parset
+)

src/geofaultnetworks/results/rock/plot.tex

+\documentclass[a4paper]{article}
+
+\usepackage{tikz}
+
+
+\begin{document}
+	\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

src/geofaultnetworks/rockfaultnetwork.cc

+#include <config.h>
+
+#include <cmath>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <string>
+
+/*
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+#include <dune/grid/uggrid/ugwrapper.hh>
+#pragma GCC diagnostic pop
+*/
+
+// dune-common includes
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/shared_ptr.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/timer.hh>
+
+// dune-istl includes
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+// dune-fufem includes
+#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
+
+// dune-grid includes
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+// dune-solver includes
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/norms/twonorm.hh>
+
+// dune-faultnetworks includes
+#include <dune/faultnetworks/utils/debugutils.hh>
+#include <dune/faultnetworks/utils/matrixreader.hh>
+#include <dune/faultnetworks/utils/vectorreader.hh>
+#include <dune/faultnetworks/utils/levelinterfacenetworkwriter.hh>
+
+#include <dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh>
+#include <dune/faultnetworks/solvers/osccgsolver.hh>
+
+#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
+#include <dune/faultnetworks/assemblers/osclocalassembler.hh>
+
+#include <dune/faultnetworks/oscrhs.hh>
+#include <dune/faultnetworks/oscdata.hh>
+#include <dune/faultnetworks/oscdatahierarchy.hh>
+
+#include <dune/faultnetworks/faultp1nodalbasis.hh>
+#include <dune/faultnetworks/interfacenetwork.hh>
+#include <dune/faultnetworks/levelinterfacenetwork.hh>
+#include <dune/faultnetworks/levelinterfacenetworkproblem.hh>
+#include <dune/faultnetworks/dgmgtransfer.hh>
+
+#include <dune/faultnetworks/faultfactories/rockfaultfactory.hh>
+
+#include <dune/solvers/solvers/loopsolver.hh>
+
+int main(int argc, char** argv) { try
+{
+    MPIHelper::instance(argc, argv);
+
+    const int dim = 2;
+    //typedef double ctype;
+
+
+    typedef typename Dune::FieldVector<double, dim> WorldVectorType;
+    typedef typename Dune::BlockVector<Dune::FieldVector<double,1>> VectorType;
+    typedef typename Dune::BitSetVector<1> BitVector;
+    typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> MatrixType;
+
+
+    #if HAVE_UG
+    typedef typename Dune::UGGrid<dim> GridType;
+    #else
+    #error No UG!
+    #endif
+
+    typedef typename GridType::LevelGridView GV;
+    typedef FaultP1NodalBasis<GV, double> DGBasis;
+    typedef typename DGBasis::LocalFiniteElement DGFE;
+
+    typedef OscLocalAssembler<GridType, DGFE, DGFE> LocalOscAssembler;
+    typedef IntersectionJumpMassAssembler<GV, DGFE, DGFE> LocalInterfaceAssembler;
+ 
+  // parse parameter file
+    ParameterTree parameterSet;
+    if (argc==2)
+        ParameterTreeParser::readINITree(argv[1], parameterSet);
+    else
+        ParameterTreeParser::readINITree("rockfaultnetwork.parset", parameterSet);
+
+    const std::string path = parameterSet.get<std::string>("path");
+    const std::string resultPath = parameterSet.get<std::string>("resultPath");
+
+    int problemCount = 0;
+    while (parameterSet.hasSub("problem" + std::to_string(problemCount))) {
+	
+    ParameterTree& problemParameters = parameterSet.sub("problem" + std::to_string(problemCount));
+
+    const std::string oscDataFile = problemParameters.get<std::string>("oscDataFile");
+
+    const int coarseResolution = problemParameters.get<int>("coarseResolution");
+    const int fineResolution = problemParameters.get<int>("fineResolution");
+    const int exactResolution = problemParameters.get<int>("exactResolution");
+
+    const double c = problemParameters.get<double>("penaltyFactor");
+
+    const int maxIterations = problemParameters.get<int>("maxIterations");
+    const double solverTolerance = problemParameters.get<double>("tol");
+    const bool nestedIteration = problemParameters.get<bool>("nestedIteration");
+
+    const int fineGridN = std::pow(2, fineResolution);
+
+    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(coarseResolution) + "_" + std::to_string(fineResolution) + "_" + std::to_string(exactResolution) + ".log");
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt
+
+
+    std::cout << std::endl;
+    std::cout << "Parameter file read successfully!" << std::endl;
+
+    typedef Dune::Matrix<Dune::FieldMatrix<double,1,1>, std::allocator<Dune::FieldMatrix<double,1,1>>> OscMatrixType;
+    OscMatrixType matrix;
+
+    // interface integral jump penalty, B:
+    OscMatrixType B(2,2);
+    B[0][0] = 1;
+    B[1][1] = 1;
+    B[0][1] = -1;
+    B[1][0] = -1;
+
+    MatrixReader<OscMatrixType> matrixReader(path + oscDataFile);
+    matrixReader.read(matrix);
+
+    const int oscGridN = matrix.N();
+    const int oscGridLevelIdx = std::max(std::round(std::log2(oscGridN)) - coarseResolution, 0.0);
+
+    if (oscGridN>fineGridN)
+        DUNE_THROW(Dune::Exception, "Provided oscData too large!");
+
+    std::cout << "OscData file read successfully!" << std::endl << std::endl;
+
+    const int fineLevelIdx = fineResolution - coarseResolution;
+    const int exactLevelIdx = exactResolution - coarseResolution;
+
+    // build multilevel fault network
+    const auto& networkParset = parameterSet.sub("network");
+    const unsigned int randomSeed = networkParset.get<unsigned int>("randomSeed");
+    const double maxAngle = networkParset.get<double>("maxAngle");
+
+    std::cout << "Creating interface network..." << std::flush;
+
+    std::srand(randomSeed);
+    RockFaultFactory<GridType> faultFactory(coarseResolution, exactLevelIdx, maxAngle);
+    const auto& interfaceNetwork = faultFactory.interfaceNetwork();
+
+    std::cout << "done!" << std::endl;
+
+    if (problemCount==0) {
+        std::vector<int> writeLevels(exactLevelIdx+1); // {0, 1, 2};
+        for (size_t i=0; i<writeLevels.size(); i++) {
+            writeLevels[i] = i;
+        }
+
+        for (size_t i=0; i<writeLevels.size(); i++) {
+            int writeLevel = writeLevels[i];
+            bool writeGrid = false; //coarseResolution+writeLevel<7;
+
+            LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
+            networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid);
+        }
+    }
+
+    const GridType& grid = faultFactory.grid();
+
+    OscData<GridType, OscMatrixType> oscData(grid, oscGridLevelIdx);
+    oscData.set(matrix);
+    OscDataHierarchy<GridType, OscMatrixType> oscDataHierarchy(oscData);
+    
+    Timer totalTimer;
+    totalTimer.reset();
+    totalTimer.start();
+
+
+    // ----------------------------------------------------------------
+    // ---              compute initial iterate                     ---
+    // ----------------------------------------------------------------
+    VectorType coarseSol;
+    std::shared_ptr<DGBasis> coarseBasis;
+
+    OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
+    L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);
+
+    std::cout << std::endl;
+    std::cout << "Computing initial iterate!" << std::endl;
+
+    if (!nestedIteration) {
+        const LevelInterfaceNetwork<GV>& coarseLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(0);
+        LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> coarseGlobalAssembler(coarseLevelInterfaceNetwork);
+        LocalOscAssembler coarseLocalAssembler(oscDataHierarchy[0]->data(), oscDataHierarchy[0]->mapper());
+
+        std::vector<std::shared_ptr<LocalInterfaceAssembler>> coarseLocalInterfaceAssemblers(coarseLevelInterfaceNetwork.size());
+        for (size_t i=0; i<coarseLocalInterfaceAssemblers.size(); i++) {
+            const int k = coarseResolution + coarseLevelInterfaceNetwork.getIntersectionsLevels().at(i);
+            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k);
+            coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
+        }
+
+        coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
+        coarseGlobalAssembler.solve();
+
+        coarseSol = coarseGlobalAssembler.getSol();
+        coarseBasis = coarseGlobalAssembler.basis();
+    } else {
+        //TODO
+        const LevelInterfaceNetwork<GV>& nLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(fineLevelIdx-1);
+        coarseBasis = std::make_shared<DGBasis>(nLevelInterfaceNetwork);
+
+        coarseSol.resize(coarseBasis->size());
+
+        std::stringstream solFilename;
+        solFilename << resultPath << "solutions/level_" << std::to_string(fineLevelIdx-1);
+
+        std::ifstream file(solFilename.str().c_str(), std::ios::in|std::ios::binary);
+        if (not(file))
+            DUNE_THROW(SolverError, "Couldn't open file " << solFilename.str() << " for reading");
+
+        Dune::MatrixVector::Generic::readBinary(file, coarseSol);
+
+        file.close();
+    }
+
+    // ----------------------------------------------------------------
+    // ---              compute exact solution                      ---
+    // ----------------------------------------------------------------
+    std::cout << std::endl;
+    std::cout << "Computing exact solution!" << std::endl;
+
+    const LevelInterfaceNetwork<GV>& exactLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(exactLevelIdx);
+    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> exactGlobalAssembler(exactLevelInterfaceNetwork);
+    LocalOscAssembler exactLocalAssembler(oscDataHierarchy[exactLevelIdx]->data(), oscDataHierarchy[exactLevelIdx]->mapper());
+
+    std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size());
+    for (size_t i=0; i<exactLocalInterfaceAssemblers.size(); i++) {
+        const int k = coarseResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(i);
+        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k);
+        exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
+    }
+
+    exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
+    exactGlobalAssembler.solve();
+
+    const VectorType& exactSol = exactGlobalAssembler.getSol();
+    std::shared_ptr<DGBasis> exactBasis = exactGlobalAssembler.basis();
+
+    EnergyNorm<MatrixType, VectorType> exactEnergyNorm(exactGlobalAssembler.matrix());
+
+    // ----------------------------------------------------------------
+    // ---                set up cg solver                          ---
+    // ----------------------------------------------------------------
+    std::cout << std::endl;
+    std::cout << "---------------------------------------------" << std::endl; 
+    std::cout << "Setting up the fine level CG solver:" << std::endl; 
+    std::cout << "---------------------------------------------" << std::endl << std::endl;
+
+    const LevelInterfaceNetwork<GV>& fineLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(fineLevelIdx);
+
+    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> fineGlobalAssembler(fineLevelInterfaceNetwork);
+    LocalOscAssembler fineLocalAssembler(oscDataHierarchy[fineLevelIdx]->data(), oscDataHierarchy[fineLevelIdx]->mapper());
+
+    std::vector<std::shared_ptr<LocalInterfaceAssembler>> fineLocalInterfaceAssemblers(fineLevelInterfaceNetwork.size());
+    for (size_t i=0; i<fineLocalInterfaceAssemblers.size(); i++) {
+        const int k = coarseResolution + fineLevelInterfaceNetwork.getIntersectionsLevels().at(i);
+        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k);
+        fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
+    }
+
+    fineGlobalAssembler.assemble(fineLocalAssembler, fineLocalInterfaceAssemblers, l2FunctionalAssembler);
+    fineGlobalAssembler.solve();
+
+    const VectorType& fineSol = fineGlobalAssembler.getSol();
+    const VectorType& fineRhs = fineGlobalAssembler.rhs();
+    std::shared_ptr<DGBasis> fineBasis = fineGlobalAssembler.basis();
+
+    EnergyNorm<MatrixType, VectorType> fineEnergyNorm(fineGlobalAssembler.matrix());
+
+    std::cout << "Setting up the preconditioner!" << std::endl;
+
+    BitVector activeLevels(fineLevelIdx+1, true);
+
+    std::vector<std::shared_ptr<LocalOscAssembler>> localAssemblers;
+    std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>> localIntersectionAssemblers;
+
+    localAssemblers.resize(activeLevels.size());
+    localIntersectionAssemblers.resize(activeLevels.size());
+
+    for (size_t i=0; i<activeLevels.size(); i++) {
+        localAssemblers[i] = std::make_shared<LocalOscAssembler>(oscDataHierarchy[i]->data(), oscDataHierarchy[i]->mapper());
+
+        const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(i);
+        std::vector<std::shared_ptr<LocalInterfaceAssembler>>& levelLocalIntersectionAssemblers = localIntersectionAssemblers[i];
+        levelLocalIntersectionAssemblers.resize(levelInterfaceNetwork.size());
+        for (size_t j=0; j<levelLocalIntersectionAssemblers.size(); j++) {
+            const int k = coarseResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j);
+            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k);
+            levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
+        }
+    }
+
+    const auto& preconditionerParset = parameterSet.sub("preconditioner");
+    using MultilevelPatchPreconditioner = MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+
+    MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerParset);
+    /*for (size_t i=0; i<preconditioner.size()-1; i++) {
+        auto& levelFaultPreconditioner = *preconditioner.levelPatchPreconditioner(i);
+
+        levelFaultPreconditioner.setPatchDepth(patchDepth);
+        levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
+    }*/
+    preconditioner.build();
+
+    std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
+    std::cout << std::endl << std::endl;
+
+    // set initial iterate
+    VectorType initialX;
+    DGMGTransfer<DGBasis> coarseFineTransfer(*coarseBasis, *fineBasis);
+    coarseFineTransfer.prolong(coarseSol, initialX);
+
+    VectorType rhs(fineRhs);
+
+    DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis);
+
+    // solve
+    OscCGSolver<MatrixType, VectorType, DGBasis>
+    solver(*fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner,
+                    maxIterations, solverTolerance, &exactEnergyNorm,
+                    Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)
+
+    /*VectorType x = initialX;
+    preconditioner.setProblem(fineGlobalAssembler.matrix(), x, rhs);
+
+    Dune::Solvers::LoopSolver<VectorType>
+    solver (preconditioner, maxIterations, solverTolerance,
+                   fineEnergyNorm,
+                   Solver::FULL,
+                   true, &fineSol); */
+
+
+    //solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec";
+
+    solver.check();
+    solver.preprocess();
+    solver.solve();
+
+    std::cout << std::endl;
+    std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl; 
+
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+
+    std::cout << "Problem " << problemCount << " done!" << std::endl;
+    problemCount++;
+    }
+} catch (Dune::Exception e) {
+
+    std::cout << e << std::endl;
+
+}
+}

src/geofaultnetworks/rockfaultnetwork.parset

+path        = ../data/
+resultPath  = ../geofaultnetworks/results/rock/
+
+[network]
+randomSeed = 5
+maxAngle = 2.0
+
+[preconditioner]
+patch = SUPPORT            # CELL , SUPPORT
+mode = ADDITIVE         # ADDITIVE , MULTIPLICATIVE
+multDirection = SYMMETRIC    # SYMMETRIC , FORWARD , BACKWARD
+patchDepth = 1
+
+###########################################
+
+[problem0]
+oscDataFile         = oscDataLaplace4.mat
+
+# level resolution in 2^(-...)
+coarseResolution = 5
+fineResolution = 8
+exactResolution = 9
+
+penaltyFactor = 1
+
+# cg solver
+nestedIteration = 0
+tol = 1e-12
+maxIterations = 8
+
+###########################################