diff --git a/dune/faultnetworks/levelinterfacenetworkproblem.hh b/dune/faultnetworks/levelinterfacenetworkproblem.hh
index ffa55382f9082822e43d5416d89ea58411771a5d..2cc53f7088394423587b422edb34e987d65a234f 100644
--- a/dune/faultnetworks/levelinterfacenetworkproblem.hh
+++ b/dune/faultnetworks/levelinterfacenetworkproblem.hh
@@ -24,7 +24,7 @@ class LevelInterfaceNetworkProblem
 
         const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
 
-        DGBasis basis_;
+        std::shared_ptr<DGBasis> basis_;
 
         MatrixType matrix_;
         VectorType rhs_;
@@ -35,7 +35,7 @@ class LevelInterfaceNetworkProblem
     public:
         LevelInterfaceNetworkProblem(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork) :
             levelInterfaceNetwork_(levelInterfaceNetwork),
-            basis_(levelInterfaceNetwork_),
+            basis_(std::make_shared<DGBasis>(levelInterfaceNetwork_)),
             requireAssemble_(true) {}
 
 
@@ -43,7 +43,7 @@ class LevelInterfaceNetworkProblem
         void assemble(LocalAssembler& localAssembler, const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers, LocalFunctionalAssembler& localFunctionalAssembler) {
             if (requireAssemble_) {
 
-                GlobalFaultAssembler<DGBasis, DGBasis> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
+                GlobalFaultAssembler<DGBasis, DGBasis> globalFaultAssembler(*basis_, *basis_, levelInterfaceNetwork_);
                 globalFaultAssembler.assembleOperator(localAssembler, localInterfaceAssemblers, matrix_);
                 globalFaultAssembler.assembleFunctional(localFunctionalAssembler, rhs_);
 
@@ -57,7 +57,7 @@ class LevelInterfaceNetworkProblem
                 // set boundary conditions
                 Dune::BitSetVector<1> boundaryDofs;
                 BoundaryPatch<GridView> boundaryPatch(gridView, true);
-                constructBoundaryDofs(boundaryPatch, basis_, boundaryDofs);
+                constructBoundaryDofs(boundaryPatch, *basis_, boundaryDofs);
 
                 for(size_t i=0; i<boundaryDofs.size(); i++) {
                     if(!boundaryDofs[i][0])
@@ -122,7 +122,7 @@ class LevelInterfaceNetworkProblem
             return x_;
         }
 
-        const DGBasis& basis() const {
+        std::shared_ptr<DGBasis> basis() const {
             if (requireAssemble_)
                 DUNE_THROW(Dune::Exception, "GlobalProblemAssembler::basis() Call assemble() before accessing basis()!");
             else
diff --git a/dune/faultnetworks/solvers/osccgsolver.cc b/dune/faultnetworks/solvers/osccgsolver.cc
index 4e52da6311ba7de32cce9d34c958d52c7215c8d3..2d6c841eaa33a69ecf30919fa34e29d2a2e1eb74 100644
--- a/dune/faultnetworks/solvers/osccgsolver.cc
+++ b/dune/faultnetworks/solvers/osccgsolver.cc
@@ -1,6 +1,8 @@
 #include <cmath>
 #include <limits>
+#include <iostream>
 #include <iomanip>
+#include <fstream>
 
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/faultnetworks/utils/debugutils.hh>
@@ -17,6 +19,21 @@ void OscCGSolver<MatrixType, VectorType, DGBasisType>::check() const
 
 }
 
+template <class MatrixType, class VectorType, class DGBasisType>
+void OscCGSolver<MatrixType, VectorType, DGBasisType>::writeSolution(const VectorType& iterate) const
+{
+    std::stringstream solFilename;
+    solFilename << this->historyBuffer_;
+
+    std::ofstream file(solFilename.str().c_str(), std::ios::out|std::ios::binary);
+    if (not(file))
+        DUNE_THROW(SolverError, "Couldn't open file " << solFilename.str() << " for writing");
+
+    Dune::MatrixVector::Generic::writeBinary(file, iterate);
+
+    file.close();
+}
+
 template <class MatrixType, class VectorType, class DGBasisType>
 void OscCGSolver<MatrixType, VectorType, DGBasisType>::solve()
 {
@@ -139,10 +156,6 @@ void OscCGSolver<MatrixType, VectorType, DGBasisType>::solve()
         p += q;                 // orthogonalization with correction
         rholast = rho;          // remember rho for recurrence
 
-        // write iteration to file, if requested
-        if (this->historyBuffer_!="") 
-            this->writeIterate(*x_, i);
-
         // Compute error
         VectorType currentDiscretizationError;
         mgTransfer_.prolong(*x_, currentDiscretizationError);
@@ -207,6 +220,10 @@ void OscCGSolver<MatrixType, VectorType, DGBasisType>::solve()
 
     }
 
+    // write last iterate to file, if requested
+    if (this->historyBuffer_!="")
+        this->writeSolution(*x_);
+
 
     if (this->verbosity_ != NumProc::QUIET) {
       std::cout << "maxTotalConvRate: " << this->maxTotalConvRate_ << ",   "
@@ -214,6 +231,8 @@ void OscCGSolver<MatrixType, VectorType, DGBasisType>::solve()
                 << i << " iterations performed\n";
         std::cout << "--------------------\n";
 
+
+
 	if (this->useRelativeError_)
             std::cout << "Norm of final discretization error: " << (discretizationError*normOfInitialDiscretizationError) << std::endl;  
         else
diff --git a/dune/faultnetworks/solvers/osccgsolver.hh b/dune/faultnetworks/solvers/osccgsolver.hh
index d91a4282a56a6e762370e6030f0bf77f2fd4ae85..1104f7a08a88aa833394aebc4562383f563f99e8 100644
--- a/dune/faultnetworks/solvers/osccgsolver.hh
+++ b/dune/faultnetworks/solvers/osccgsolver.hh
@@ -48,6 +48,8 @@ public:
      */
     virtual void check() const;
 
+    void writeSolution(const VectorType& iterate) const;
+
     const DGBasisType& fineBasis_;
 
     const MatrixType* matrix_;
diff --git a/dune/faultnetworks/utils/levelinterfacenetworkwriter.hh b/dune/faultnetworks/utils/levelinterfacenetworkwriter.hh
index b4f43e15c5417d324c286ba05535f9c2858b24b6..be074b5819b3adf8a04e170e26efbb00cfab242c 100644
--- a/dune/faultnetworks/utils/levelinterfacenetworkwriter.hh
+++ b/dune/faultnetworks/utils/levelinterfacenetworkwriter.hh
@@ -3,6 +3,7 @@
 
 #include <string>
 #include <dune/common/fvector.hh>
+#include <dune/grid/common/mcmgmapper.hh>
 #include <dune/fufem/boundaryiterator.hh>
 //#include <dune/faultnetworks/levelinterfacenetwork.hh>
 
@@ -10,8 +11,18 @@ class LevelInterfaceNetworkWriter {
     private:
         std::string filePath_;
 
+        //! Parameter for mapper class
+        template<int dim>
+        struct FaceMapperLayout
+        {
+            bool contains (Dune::GeometryType gt)
+            {
+                return gt.dim() == dim-1;
+            }
+        };
+
         template <typename Intersection>
-        void writeIntersection(const Intersection& intersection, std::ofstream& file, const int colorID) {
+        void writeIntersection(const Intersection& intersection, std::ofstream& file, const std::string colorTag) {
             typedef typename Dune::FieldVector<typename Intersection::ctype, Intersection::dimensionworld > GlobalCoordinates;
 
             const auto& geometry = intersection.geometry();
@@ -26,27 +37,7 @@ class LevelInterfaceNetworkWriter {
 
             for (size_t i=0; i<vertices.size(); i++) {
                 for (size_t j=i+1; j<vertices.size(); j++) {
-
-                    switch (colorID) {
-                        case 0:
-                            file << "\\draw[black] ";
-                            break;
-                        case 1:
-                            file << "\\draw[MATLABblue] ";
-                            break;
-                        case 2:
-                            file << "\\draw[MATLABgreen] ";
-                            break;
-                        case 3:
-                            file << "\\draw[MATLABred] ";
-                            break;
-                        case 4:
-                            file << "\\draw[MATLAByellow] ";
-                            break;
-                        default:
-                            file << "\\draw[black] ";
-                    }
-
+                    file << "\\draw[" << colorTag << "] ";
                     writeVertex(vertices[i], file);
                     file << " -- ";
                     writeVertex(vertices[j], file);
@@ -77,7 +68,7 @@ class LevelInterfaceNetworkWriter {
 
         // works only in 1D and 2D
         template <typename LevelInterfaceNetwork>
-        void write(const LevelInterfaceNetwork& levelInterfaceNetwork, bool highlightNewLevelInterfaces = false) {
+        void write(const LevelInterfaceNetwork& levelInterfaceNetwork, bool writeGrid = false) {
 			typedef typename LevelInterfaceNetwork::GridType::LevelGridView GridView;
             typedef typename LevelInterfaceNetwork::Intersection Intersection;
 
@@ -90,26 +81,69 @@ class LevelInterfaceNetworkWriter {
                 file << "\\definecolor{MATLABgreen}{RGB}{60,140,40}\n";
                 file << "\\definecolor{MATLABorange}{RGB}{212,83,25}\n";
                 file << "\\definecolor{MATLAByellow}{RGB}{237,177,32}\n";
+                file << "\\definecolor{lightGray}{RGB}{211,211,211}\n";
+                file << "\n";
+
+                file << "\\colorlet{boundary}{black}\n";
+                file << "\\colorlet{grid}{lightGray}\n";
+                for (size_t i=0; i<=levelInterfaceNetwork.level(); i++) {
+
+                    file << "\\colorlet{level" << i << "}{black}\n";
+                }
+                file << "\n";
 
                 file << "\\begin{tikzpicture}[scale=\\scale]\n";
 
                 const GridView& gridView = levelInterfaceNetwork.levelGridView();
 
-                // domain boundary
-                BoundaryIterator<GridView> bIt(gridView, BoundaryIterator<GridView>::begin);
-                BoundaryIterator<GridView> bEnd(gridView, BoundaryIterator<GridView>::end);
-				for(; bIt!=bEnd; ++bIt) {
-                    writeIntersection(*bIt, file, 0);
-				}
-
-				// fault intersections
-                const std::vector<Intersection>& intersections = levelInterfaceNetwork.getIntersections();
-                const std::vector<int>& intersectionsLevels = levelInterfaceNetwork.getIntersectionsLevels();
-                for (size_t i=0; i<intersections.size(); i++) {
-                    if (highlightNewLevelInterfaces) {
-                        writeIntersection(intersections[i], file, 3*(intersectionsLevels[i]==levelInterfaceNetwork.level()));
-                    } else {
-                        writeIntersection(intersections[i], file, intersectionsLevels[i]);
+                std::string colorTag;
+
+                if (writeGrid) {
+                    typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView, FaceMapperLayout > FaceMapper;
+                    FaceMapper intersectionMapper(gridView);
+                    std::vector<bool> intersectionHandled(intersectionMapper.size(),false);
+
+                    // boundary and grid intersections
+                    for (const auto& elem:elements(gridView)) {
+                        for (const auto& isect:intersections(gridView, elem)) {
+
+                            if (intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)])
+                                continue;
+
+                            intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)] = true;
+
+                            if (isect.boundary()) {
+                                colorTag = "boundary, very thin";
+                            } else if (!levelInterfaceNetwork.isInterfaceIntersection(isect)) {
+                                colorTag = "grid, very thin";
+                            }
+                            writeIntersection(isect, file, colorTag);
+                        }
+                    }
+
+                    // fault intersections
+                    const std::vector<Intersection>& intersections = levelInterfaceNetwork.getIntersections();
+                    const std::vector<int>& intersectionsLevels = levelInterfaceNetwork.getIntersectionsLevels();
+                    for (size_t i=0; i<intersections.size(); i++) {
+                        colorTag = "level" + std::to_string(intersectionsLevels[i]);
+                        writeIntersection(intersections[i], file, colorTag);
+                    }
+
+                } else {
+                    // domain boundary
+                    BoundaryIterator<GridView> bIt(gridView, BoundaryIterator<GridView>::begin);
+                    BoundaryIterator<GridView> bEnd(gridView, BoundaryIterator<GridView>::end);
+                    for(; bIt!=bEnd; ++bIt) {
+                        colorTag = "boundary, very thin";
+                        writeIntersection(*bIt, file, colorTag);
+                    }
+
+                    // fault intersections
+                    const std::vector<Intersection>& intersections = levelInterfaceNetwork.getIntersections();
+                    const std::vector<int>& intersectionsLevels = levelInterfaceNetwork.getIntersectionsLevels();
+                    for (size_t i=0; i<intersections.size(); i++) {
+                        colorTag = "level" + std::to_string(intersectionsLevels[i]);
+                        writeIntersection(intersections[i], file, colorTag);
                     }
                 }
 
diff --git a/src/cantorfaultnetworks/cantorconvergence.cc b/src/cantorfaultnetworks/cantorconvergence.cc
index 271b9fbdedd418317da988b66a473c40c3a645e6..a6769fd56c7be64bbd9b274bd818442ea16d09fb 100644
--- a/src/cantorfaultnetworks/cantorconvergence.cc
+++ b/src/cantorfaultnetworks/cantorconvergence.cc
@@ -43,7 +43,6 @@
 #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>
@@ -153,6 +152,14 @@ int main(int argc, char** argv) { try
     CantorConvergenceFaultFactory<GridType> faultFactory(minResolution, exactLevelIdx, minCantorResolution, maxCantorLevels);
     const std::vector<std::shared_ptr<InterfaceNetwork<GridType>>>& interfaceNetworks = faultFactory.interfaceNetworks();
 
+    /*
+    const InterfaceNetwork<GridType>& interfaceNetwork = *interfaceNetworks[interfaceNetworks.size()-1];
+       for (size_t i=0; i<interfaceNetwork.size(); i++) {
+           LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
+           networkWriter.write(interfaceNetwork.levelInterfaceNetwork(i), false);
+       }
+    */
+
     std::cout << "Interface networks built successfully!" << std::endl << std::endl;
 
     const GridType& grid = faultFactory.grid();
diff --git a/src/cantorfaultnetworks/cantorconvergence.parset b/src/cantorfaultnetworks/cantorconvergence.parset
index 7944ff84a6630c9942975f76eb22580205ab83aa..a7a874f6247089763a816ce631a63d64adedb0bc 100644
--- a/src/cantorfaultnetworks/cantorconvergence.parset
+++ b/src/cantorfaultnetworks/cantorconvergence.parset
@@ -1,5 +1,5 @@
 path        = ../data/
-resultPath  = ../results/cantorfaultnetworks/cantorconvergence/
+resultPath  = ../cantorfaultnetworks/results/cantorconvergence/
 
 
 [problem0]
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.cc b/src/cantorfaultnetworks/cantorfaultnetwork.cc
index af5dd15693e161a8a7032d1090582865eea07091..eaf5a8c89b3ec67ec1e1c8bf356bdbd9930f591f 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.cc
@@ -80,7 +80,6 @@ int main(int argc, char** argv) { try
     #error No UG!
     #endif
 
-
     typedef typename GridType::LevelGridView GV;
     typedef FaultP1NodalBasis<GV, double> DGBasis;
     typedef typename DGBasis::LocalFiniteElement DGFE;
@@ -115,6 +114,7 @@ int main(int argc, char** argv) { try
 
     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);
 
@@ -160,10 +160,13 @@ int main(int argc, char** argv) { try
     CantorFaultFactory<GridType> faultFactory(levelToCantorLevel, coarseResolution, exactLevelIdx, maxCantorLevel);
     const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
 
-    if (problemCount==0) {
-        int writeLevel = fineLevelIdx;
-        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork.tikz");
-        networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), false);
+    std::vector<int> writeLevels {0, 1, fineLevelIdx};
+    for (size_t i=0; i<writeLevels.size(); i++) {
+        int writeLevel = writeLevels[i];
+        bool writeGrid = !(writeLevel==fineLevelIdx);
+
+        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
+        networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid);
     }
 
     const GridType& grid = faultFactory.grid();
@@ -180,28 +183,49 @@ int main(int argc, char** argv) { try
     // ----------------------------------------------------------------
     // ---              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;
 
-    const LevelInterfaceNetwork<GV>& coarseLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(0);
-    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> coarseGlobalAssembler(coarseLevelInterfaceNetwork);
-    LocalOscAssembler coarseLocalAssembler(oscDataHierarchy[0]->data(), oscDataHierarchy[0]->mapper());
+    if (!nestedIteration) {
+        const LevelInterfaceNetwork<GV>& coarseLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(0);
+        LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> coarseGlobalAssembler(coarseLevelInterfaceNetwork);
+        LocalOscAssembler coarseLocalAssembler(oscDataHierarchy[0]->data(), oscDataHierarchy[0]->mapper());
 
-    OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
-    L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);
+        std::vector<std::shared_ptr<LocalInterfaceAssembler>> coarseLocalInterfaceAssemblers(coarseLevelInterfaceNetwork.size());
+        for (size_t i=0; i<coarseLocalInterfaceAssemblers.size(); i++) {
+            const int k = minCantorResolution + coarseLevelInterfaceNetwork.getIntersectionsLevels().at(i);
+            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
+            coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
+        }
 
-    std::vector<std::shared_ptr<LocalInterfaceAssembler>> coarseLocalInterfaceAssemblers(coarseLevelInterfaceNetwork.size());
-    for (size_t i=0; i<coarseLocalInterfaceAssemblers.size(); i++) {
-        const int k = minCantorResolution + coarseLevelInterfaceNetwork.getIntersectionsLevels().at(i);
-        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
-        coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
-    }
+        coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
+        coarseGlobalAssembler.solve();
+
+        coarseSol = coarseGlobalAssembler.getSol();
+        coarseBasis = coarseGlobalAssembler.basis();
+    } else {
+        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(fineResolution-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");
 
-    coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
-    coarseGlobalAssembler.solve();
+        Dune::MatrixVector::Generic::readBinary(file, coarseSol);
 
-    const VectorType& coarseSol = coarseGlobalAssembler.getSol();
-    const DGBasis& coarseBasis = coarseGlobalAssembler.basis();
+        file.close();
+    }
 
     // ----------------------------------------------------------------
     // ---              compute exact solution                      ---
@@ -224,7 +248,7 @@ int main(int argc, char** argv) { try
     exactGlobalAssembler.solve();
 
     const VectorType& exactSol = exactGlobalAssembler.getSol();
-    const DGBasis& exactBasis = exactGlobalAssembler.basis();
+    std::shared_ptr<DGBasis> exactBasis = exactGlobalAssembler.basis();
 
     EnergyNorm<MatrixType, VectorType> exactEnergyNorm(exactGlobalAssembler.matrix());
 
@@ -253,13 +277,10 @@ int main(int argc, char** argv) { try
 
     const VectorType& fineSol = fineGlobalAssembler.getSol();
     const VectorType& fineRhs = fineGlobalAssembler.rhs();
-    const DGBasis& fineBasis = fineGlobalAssembler.basis();
+    std::shared_ptr<DGBasis> fineBasis = fineGlobalAssembler.basis();
 
     EnergyNorm<MatrixType, VectorType> fineEnergyNorm(fineGlobalAssembler.matrix());
 
-    // multilevel transfer
-    DGMGTransfer<DGBasis> coarseFineTransfer(coarseBasis, fineBasis);
-
     std::cout << "Setting up the preconditioner!" << std::endl;
 
     BitVector activeLevels(fineLevelIdx+1, true);
@@ -299,18 +320,24 @@ int main(int argc, char** argv) { try
     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);
+    DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis);
 
     // solve
     OscCGSolver<MatrixType, VectorType, DGBasis>
-    solver(fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner,
+    solver(*fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner,
                     maxIterations, solverTolerance, &exactEnergyNorm,
                     Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)
+
+
+    solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineResolution);
+
     solver.check();
     solver.preprocess();
     solver.solve();
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.parset b/src/cantorfaultnetworks/cantorfaultnetwork.parset
index d17010a5d687f98d331b3ec68feb9f49e6cf5dbc..0d8411bb631bf9e36af423302f000bd1d9312dc8 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.parset
@@ -1,7 +1,8 @@
 path        = ../data/
-resultPath  = ../results/cantorfaultnetworks/
+resultPath  = ../cantorfaultnetworks/results/
 
 
+###########################################
 
 [problem0]
 oscDataFile         = oscDataLaplace32.mat
@@ -16,6 +17,7 @@ penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
@@ -34,6 +36,7 @@ penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
@@ -52,6 +55,7 @@ penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
@@ -70,6 +74,7 @@ penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
@@ -89,6 +94,7 @@ penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
diff --git a/src/results/.gitkeep b/src/cantorfaultnetworks/results/cantorconvergence/.gitkeep
similarity index 100%
rename from src/results/.gitkeep
rename to src/cantorfaultnetworks/results/cantorconvergence/.gitkeep
diff --git a/src/data/oscDataLaplace4.mat b/src/data/oscDataLaplace4.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b40a5e81fb8b25ffa4c89a31586d7fb63ddc96a0
--- /dev/null
+++ b/src/data/oscDataLaplace4.mat
@@ -0,0 +1,5 @@
+4 4 
+1 1 1 1
+1 1 1 1
+1 1 1 1
+1 1 1 1
diff --git a/src/geofaultnetworks/geofaultnetwork.cc b/src/geofaultnetworks/geofaultnetwork.cc
index 5dadf92fbd4f0f775c930f7fd5da4c485022ac7c..99b5ace89c01f984d3f422cd3dcdf4c45e7e17d9 100644
--- a/src/geofaultnetworks/geofaultnetwork.cc
+++ b/src/geofaultnetworks/geofaultnetwork.cc
@@ -63,7 +63,6 @@
 
 #include <dune/faultnetworks/faultfactories/geofaultfactory.hh>
 
-
 int main(int argc, char** argv) { try
 {
     MPIHelper::instance(argc, argv);
@@ -119,6 +118,7 @@ int main(int argc, char** argv) { try
 
     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);
 
@@ -180,6 +180,9 @@ int main(int argc, char** argv) { try
             faultToLevel[i] = j;
             i++;
         }
+
+        faultToLevel[i] = fineLevelIdx;
+        i++;
     }
 
     fineLevelIdx = fineResolution - coarseResolution;
@@ -191,12 +194,16 @@ int main(int argc, char** argv) { try
     const GridType& grid = faultFactory.grid();
     const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
 
-    if (problemCount==0) {
-        int writeLevel = fineLevelIdx;
-        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork.tikz");
-        networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), false);
+    std::vector<int> writeLevels {0, 1, fineLevelIdx};
+    for (size_t i=0; i<writeLevels.size(); i++) {
+        int writeLevel = writeLevels[i];
+        bool writeGrid = !(writeLevel==fineLevelIdx);
+
+        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
+        networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid);
     }
 
+
     OscData<GridType, OscMatrixType> oscData(grid, oscGridLevelIdx);
     oscData.set(matrix);
     OscDataHierarchy<GridType, OscMatrixType> oscDataHierarchy(oscData);
@@ -209,28 +216,50 @@ int main(int argc, char** argv) { try
     // ----------------------------------------------------------------
     // ---              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;
 
-    const LevelInterfaceNetwork<GV>& coarseLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(0);
-    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> coarseGlobalAssembler(coarseLevelInterfaceNetwork);
-    LocalOscAssembler coarseLocalAssembler(oscDataHierarchy[0]->data(), oscDataHierarchy[0]->mapper());
+    if (!nestedIteration) {
+        const LevelInterfaceNetwork<GV>& coarseLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(0);
+        LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> coarseGlobalAssembler(coarseLevelInterfaceNetwork);
+        LocalOscAssembler coarseLocalAssembler(oscDataHierarchy[0]->data(), oscDataHierarchy[0]->mapper());
 
-    OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
-    L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);
+        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) - 1.0);
+            coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
+        }
 
-    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-1);
-        coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
-    }
+        coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
+        coarseGlobalAssembler.solve();
+
+        coarseSol = coarseGlobalAssembler.getSol();
+        coarseBasis = coarseGlobalAssembler.basis();
+    } else {
+        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(fineResolution-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");
 
-    coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
-    coarseGlobalAssembler.solve();
+        Dune::MatrixVector::Generic::readBinary(file, coarseSol);
+
+        file.close();
+    }
 
-    const VectorType& coarseSol = coarseGlobalAssembler.getSol();
-    const DGBasis& coarseBasis = coarseGlobalAssembler.basis();
 
     // ----------------------------------------------------------------
     // ---              compute exact solution                      ---
@@ -245,7 +274,7 @@ int main(int argc, char** argv) { try
     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-1);
+        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * (std::pow(2.0, k) - 1.0);
         exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
     }
 
@@ -253,7 +282,7 @@ int main(int argc, char** argv) { try
     exactGlobalAssembler.solve();
 
     const VectorType& exactSol = exactGlobalAssembler.getSol();
-    const DGBasis& exactBasis = exactGlobalAssembler.basis();
+    std::shared_ptr<DGBasis> exactBasis = exactGlobalAssembler.basis();
 
     EnergyNorm<MatrixType, VectorType> exactEnergyNorm(exactGlobalAssembler.matrix());
 
@@ -273,7 +302,7 @@ int main(int argc, char** argv) { try
     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-1);
+        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * (std::pow(2.0, k) - 1.0);
         fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
     }
 
@@ -282,13 +311,10 @@ int main(int argc, char** argv) { try
 
     const VectorType& fineSol = fineGlobalAssembler.getSol();
     const VectorType& fineRhs = fineGlobalAssembler.rhs();
-    const DGBasis& fineBasis = fineGlobalAssembler.basis();
+    std::shared_ptr<DGBasis> fineBasis = fineGlobalAssembler.basis();
 
     EnergyNorm<MatrixType, VectorType> fineEnergyNorm(fineGlobalAssembler.matrix());
 
-    // multilevel transfer
-    DGMGTransfer<DGBasis> coarseFineTransfer(coarseBasis, fineBasis);
-
     std::cout << "Setting up the preconditioner!" << std::endl;
 
     BitVector activeLevels(fineLevelIdx+1, true);
@@ -307,7 +333,7 @@ int main(int argc, char** argv) { try
         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-1);
+            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * (std::pow(2.0, k) - 1.0);
             levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
         }
     }
@@ -328,18 +354,23 @@ int main(int argc, char** argv) { try
     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);
+    DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis);
 
     // solve
     OscCGSolver<MatrixType, VectorType, DGBasis>
-    solver(fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner,
+    solver(*fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner,
                     maxIterations, solverTolerance, &exactEnergyNorm,
                     Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)
+
+    solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineResolution);
+
     solver.check();
     solver.preprocess();
     solver.solve();
diff --git a/src/geofaultnetworks/geofaultnetwork.parset b/src/geofaultnetworks/geofaultnetwork.parset
index 3bb72548ad3d5ef3c6a5254a82caa0e458600e59..9b6d5c511e1538f8c7cf9715923ec42bc430a02a 100644
--- a/src/geofaultnetworks/geofaultnetwork.parset
+++ b/src/geofaultnetworks/geofaultnetwork.parset
@@ -1,5 +1,5 @@
 path        = ../data/
-resultPath  = ../results/geofaultnetworks/
+resultPath  = ../geofaultnetworks/results/
 
 
 [problem0]
@@ -9,12 +9,13 @@ oscDataFile         = oscDataLaplace32.mat
 coarseResolution = 4
 fineResolution = 5
 exactResolution = 10
-faultCount = 30
+faultCount = 33
 
 penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
@@ -27,12 +28,13 @@ oscDataFile         = oscDataLaplace32.mat
 coarseResolution = 4
 fineResolution = 6
 exactResolution = 10
-faultCount = 30
+faultCount = 33
 
 penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
@@ -45,12 +47,13 @@ oscDataFile         = oscDataLaplace32.mat
 coarseResolution = 4
 fineResolution = 7
 exactResolution = 10
-faultCount = 30
+faultCount = 33
 
 penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
@@ -63,12 +66,13 @@ oscDataFile         = oscDataLaplace32.mat
 coarseResolution = 4
 fineResolution = 8
 exactResolution = 10
-faultCount = 30
+faultCount = 33
 
 penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
@@ -81,12 +85,13 @@ oscDataFile         = oscDataLaplace32.mat
 coarseResolution = 4
 fineResolution = 9
 exactResolution = 10
-faultCount = 30
+faultCount = 33
 
 penaltyFactor = 1
 patchDepth = 1
 
 # cg solver
+nestedIteration = 0
 tol = 1e-12
 maxIterations = 8
 
diff --git a/src/results/cantorfaultnetworks/.gitkeep b/src/results/cantorfaultnetworks/.gitkeep
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/src/results/cantorfaultnetworks/cantorconvergence/.gitkeep b/src/results/cantorfaultnetworks/cantorconvergence/.gitkeep
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/src/results/cantorfaultnetworks/plot.tex b/src/results/cantorfaultnetworks/plot.tex
deleted file mode 100644
index 5f0576d36dfcf10200cfa18a40fe4f7450244d54..0000000000000000000000000000000000000000
--- a/src/results/cantorfaultnetworks/plot.tex
+++ /dev/null
@@ -1,9 +0,0 @@
-\documentclass[tikz]{standalone}
-
-\usepackage{tikz}
-
-\begin{document}
-
-\def\scale{3}
-\input{levelinterfacenetwork_1.tikz}
-\end{document}
\ No newline at end of file
diff --git a/src/results/geofaultnetworks/.gitkeep b/src/results/geofaultnetworks/.gitkeep
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/src/results/geofaultnetworks/plot.tex b/src/results/geofaultnetworks/plot.tex
deleted file mode 100644
index b9761d1b2e3028c434b385d429265eb0356586b2..0000000000000000000000000000000000000000
--- a/src/results/geofaultnetworks/plot.tex
+++ /dev/null
@@ -1,9 +0,0 @@
-\documentclass[tikz]{standalone}
-
-\usepackage{tikz}
-
-\begin{document}
-
-\def\scale{3}
-\input{levelinterfacenetwork_5.tikz}
-\end{document}
\ No newline at end of file