.

src/CMakeLists.txt

 add_subdirectory("cantorfaultnetworks")
 add_subdirectory("geofaultnetworks")
-add_subdirectory("uniformfaultnetworks")
 
 #add_executable("multilevellaplace" multilevellaplace.cc)
 #target_link_dune_default_libraries("multilevellaplace")

src/cantorfaultnetworks/CMakeLists.txt

@@ -4,9 +4,6 @@ target_link_dune_default_libraries("cantorfaultnetwork")
 add_executable("cantorconvergence" cantorconvergence.cc)
 target_link_dune_default_libraries("cantorconvergence")
 
-add_executable("cantorconvergence2" cantorconvergence2.cc)
-target_link_dune_default_libraries("cantorconvergence2")
-
 add_custom_target(cantorfaultnetworks_srcs SOURCES
 cantorfaultnetwork.parset
 cantorconvergence.parset

src/cantorfaultnetworks/cantorconvergence.cc

@@ -58,7 +58,7 @@
 #include <dune/faultnetworks/levelinterfacenetworkproblem.hh>
 #include <dune/faultnetworks/dgmgtransfer.hh>
 
-#include <dune/faultnetworks/faultfactories/cantorfaultfactory.hh>
+#include <dune/faultnetworks/faultfactories/cantorconvergencefaultfactory.hh>
 
 int main(int argc, char** argv) { try
 {
@@ -86,7 +86,7 @@ int main(int argc, char** argv) { try
 
     typedef OscLocalAssembler<GridType, DGFE, DGFE> LocalOscAssembler;
     typedef IntersectionJumpMassAssembler<GV, DGFE, DGFE> LocalInterfaceAssembler;
- 
+
   // parse parameter file
     ParameterTree parameterSet;
     if (argc==2)
@@ -99,21 +99,22 @@ int main(int argc, char** argv) { try
 
     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 minResolution = problemParameters.get<int>("minResolution");
     const int maxResolution = problemParameters.get<int>("maxResolution");
+    const int exactResolution = problemParameters.get<int>("exactResolution");
     const int minCantorResolution = problemParameters.get<int>("minCantorResolution");
 
-    const double penaltyFactor = problemParameters.get<double>("penaltyFactor");
+    const double c = problemParameters.get<double>("penaltyFactor");
 
     const int maxGridN = std::pow(2, maxResolution);
 
 
-    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minResolution) + "_" + std::to_string(maxResolution) + "_" + std::to_string(minCantorResolution) +  "_" + std::to_string((int) penaltyFactor) + ".log");
+    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minResolution) + "_" + std::to_string(maxResolution) + "_" + std::to_string(exactResolution) + "_" + std::to_string(minCantorResolution) +  "_" + std::to_string((int) c) + ".log");
     std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
     std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt
 
@@ -143,29 +144,30 @@ int main(int argc, char** argv) { try
     std::cout << "OscData file read successfully!" << std::endl << std::endl;
 
     const int maxLevelIdx = maxResolution - minResolution;
-
-    size_t maxCantorLevel = maxLevelIdx+minCantorResolution;
-    std::map<size_t, size_t> levelToCantorLevel;
-    for (int i=0; i<=maxLevelIdx; i++) {
-        levelToCantorLevel[i] = i+minCantorResolution;
+    const int exactLevelIdx = exactResolution - minResolution;
+    std::vector<size_t> maxCantorLevels(exactLevelIdx+1);
+    for (int i=0; i<=exactLevelIdx; i++) {
+        maxCantorLevels[i] = i+minCantorResolution;
     }
 
     // build multilevel cantor fault network
-    CantorFaultFactory<GridType> faultFactory(levelToCantorLevel, minResolution, maxLevelIdx, maxCantorLevel);
-    const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
-    for (size_t i=0; i<interfaceNetwork.size(); i++) {
+    CantorConvergenceFaultFactory<GridType> faultFactory(minResolution, exactLevelIdx, minCantorResolution, maxCantorLevels);
+    const std::vector<std::shared_ptr<InterfaceNetwork<GridType>>>& interfaceNetworks = faultFactory.interfaceNetworks();
+
+    /*
+    for (size_t i=0; i<interfaceNetworks[1]->size(); i++) {
         LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
-        networkWriter.write(interfaceNetwork.levelInterfaceNetwork(i), true);
-    }
+        networkWriter.write(interfaceNetworks[1]->levelInterfaceNetwork(i), true);
+    }*/
 
-    std::cout << "Interface network built successfully!" << std::endl << std::endl;
+    std::cout << "Interface networks built successfully!" << std::endl << std::endl;
 
     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();
@@ -178,41 +180,48 @@ int main(int argc, char** argv) { try
     std::cout << "Computing direct solutions!" << std::endl;
 
     std::vector<VectorType> solutions;
-    solutions.resize(maxLevelIdx+1);
-    std::vector<double> discErrors;
-    discErrors.resize(maxLevelIdx+1);
+    solutions.resize(exactLevelIdx+1);
+
+    std::vector<double> approxErrors;
+    approxErrors.resize(maxLevelIdx+1);
+
+    double discError = NAN;
 
     OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
     L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);
 
-    const LevelInterfaceNetwork<GV>& maxLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(maxLevelIdx);
-    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> maxGlobalAssembler(maxLevelInterfaceNetwork);
-    LocalOscAssembler maxLocalAssembler(oscDataHierarchy[maxLevelIdx]->data(), oscDataHierarchy[maxLevelIdx]->mapper());
+    const LevelInterfaceNetwork<GV>& exactLevelInterfaceNetwork = interfaceNetworks.back()->levelInterfaceNetwork(exactLevelIdx);
+    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> exactGlobalAssembler(exactLevelInterfaceNetwork);
+    LocalOscAssembler exactLocalAssembler(oscDataHierarchy[exactLevelIdx]->data(), oscDataHierarchy[exactLevelIdx]->mapper());
 
-    std::vector<std::shared_ptr<LocalInterfaceAssembler>> maxLocalInterfaceAssemblers(maxLevelInterfaceNetwork.size());
-    for (size_t j=0; j<maxLocalInterfaceAssemblers.size(); j++) {
+    std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size());
+    for (size_t j=0; j<exactLocalInterfaceAssemblers.size(); j++) {
         //const double pen = penaltyFactor;
-        const double pen = penaltyFactor*std::pow(2.0, minResolution + maxLevelInterfaceNetwork.getIntersectionsLevels().at(j));
-        maxLocalInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
+        const int k = minCantorResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(j);
+        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
+        exactLocalInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
     }
 
-    maxGlobalAssembler.assemble(maxLocalAssembler, maxLocalInterfaceAssemblers, l2FunctionalAssembler);
-    maxGlobalAssembler.solve();
+    exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
+    exactGlobalAssembler.solve();
+
+    solutions[exactLevelIdx] = exactGlobalAssembler.getSol();
+    approxErrors[maxLevelIdx] = 0;
 
-    solutions[maxLevelIdx] = maxGlobalAssembler.getSol();
-    discErrors[maxLevelIdx] = 0;
+    const DGBasis& exactBasis = exactGlobalAssembler.basis();
 
-    const DGBasis& maxBasis = maxGlobalAssembler.basis();
+    EnergyNorm<MatrixType, VectorType> energyNorm(exactGlobalAssembler.matrix());
 
-    for (int i=maxLevelIdx-1; i>=0; i--) {
-        const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(i);
+    for (int i=maxLevelIdx; i>=0; i--) {
+        const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetworks[i]->levelInterfaceNetwork(maxLevelIdx);
         LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> globalAssembler(levelInterfaceNetwork);
-        LocalOscAssembler localAssembler(oscDataHierarchy[i]->data(), oscDataHierarchy[i]->mapper());
+        LocalOscAssembler localAssembler(oscDataHierarchy[maxLevelIdx]->data(), oscDataHierarchy[maxLevelIdx]->mapper());
 
         std::vector<std::shared_ptr<LocalInterfaceAssembler>> localInterfaceAssemblers(levelInterfaceNetwork.size());
         for (size_t j=0; j<localInterfaceAssemblers.size(); j++) {
             //const double pen = penaltyFactor;
-            const double pen = penaltyFactor*std::pow(2.0, minResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j));
+            const int k = minCantorResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j);
+            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
             localInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
         }
 
@@ -222,24 +231,36 @@ int main(int argc, char** argv) { try
         const VectorType& sol = globalAssembler.getSol();
         const DGBasis& basis = globalAssembler.basis();
 
-        DGMGTransfer<DGBasis> levelTransfer(basis, maxBasis);
+        DGMGTransfer<DGBasis> levelTransfer(basis, exactBasis);
         levelTransfer.prolong(sol, solutions[i]);
 
-        VectorType discError(solutions[maxLevelIdx]);
-        discError -= solutions[i];
-        discErrors[i] = discError.two_norm();
+        if (i==maxLevelIdx) {
+            VectorType discErrorV(solutions[exactLevelIdx]);
+            discErrorV -= solutions[maxLevelIdx];
+            //approxErrors[i] = approxError.two_norm();
+            discError= energyNorm.operator ()(discErrorV);
+        } else {
+            VectorType approxError(solutions[maxLevelIdx]);
+            approxError -= solutions[i];
+            //approxErrors[i] = approxError.two_norm();
+            approxErrors[i] = energyNorm.operator ()(approxError);
+        }
     }
 
     // print results
     std::cout << std::endl << std::endl;
-    std::cout << "level  discError"<< std::endl;
-    std::cout << "----------------" << std::endl;
-    for (size_t i=0; i<discErrors.size(); i++) {
-        std::cout << i << "      " << discErrors[i] << std::endl;
+    std::cout << "cantorLevel  h     approxError"<< std::endl;
+    std::cout << "------------------------------" << std::endl;
+    for (size_t i=0; i<approxErrors.size(); i++) {
+        std::cout << (minCantorResolution+i) << "            2^-" << maxResolution << "  " << approxErrors[i] << std::endl;
     }
 
+    std::cout << std::endl << std::endl;
+    std::cout << "exactResolution: 2^-" << exactResolution << std::endl;
+    std::cout << "discError: " << discError << std::endl;
+
     std::cout << std::endl;
-    std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl; 
+    std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl;
 
     std::cout.rdbuf(coutbuf); //reset to standard output again
 
@@ -247,6 +268,8 @@ int main(int argc, char** argv) { try
     problemCount++;
     }
 } catch (Dune::Exception e) {
+
     std::cout << e << std::endl;
+
 }
 }

src/cantorfaultnetworks/cantorconvergence2.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/cantorconvergencefaultfactory.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::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("/home/mi/podlesny/software/dune/faultnetworks/src/cantorfaultnetworks/cantorconvergence.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 minResolution = problemParameters.get<int>("minResolution");
-    const int maxResolution = problemParameters.get<int>("maxResolution");
-    const int exactResolution = problemParameters.get<int>("exactResolution");
-    const int minCantorResolution = problemParameters.get<int>("minCantorResolution");
-
-    const double c = problemParameters.get<double>("penaltyFactor");
-
-    const int maxGridN = std::pow(2, maxResolution);
-
-
-    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minResolution) + "_" + std::to_string(maxResolution) + "_" + std::to_string(exactResolution) + "_" + std::to_string(minCantorResolution) +  "_" + std::to_string((int) c) + ".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::round(std::log2(oscGridN)) - minResolution;
-
-    if (oscGridN>maxGridN)
-        DUNE_THROW(Dune::Exception, "Provided oscData too large!");
-
-    std::cout << "OscData file read successfully!" << std::endl << std::endl;
-
-    const int maxLevelIdx = maxResolution - minResolution;
-    const int exactLevelIdx = exactResolution - minResolution;
-    std::vector<size_t> maxCantorLevels(exactLevelIdx+1);
-    for (int i=0; i<=exactLevelIdx; i++) {
-        maxCantorLevels[i] = i+minCantorResolution;
-    }
-
-    // build multilevel cantor fault network
-    CantorConvergenceFaultFactory<GridType> faultFactory(minResolution, exactLevelIdx, minCantorResolution, maxCantorLevels);
-    const std::vector<std::shared_ptr<InterfaceNetwork<GridType>>>& interfaceNetworks = faultFactory.interfaceNetworks();
-
-    /*
-    for (size_t i=0; i<interfaceNetworks[1]->size(); i++) {
-        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
-        networkWriter.write(interfaceNetworks[1]->levelInterfaceNetwork(i), true);
-    }*/
-
-    std::cout << "Interface networks built successfully!" << std::endl << std::endl;
-
-    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 direct solutions                    ---
-    // ----------------------------------------------------------------
-    std::cout << std::endl;
-    std::cout << "Computing direct solutions!" << std::endl;
-
-    std::vector<VectorType> solutions;
-    solutions.resize(exactLevelIdx+1);
-
-    std::vector<double> approxErrors;
-    approxErrors.resize(maxLevelIdx+1);
-
-    double discError = NAN;
-
-    OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
-    L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);
-
-    const LevelInterfaceNetwork<GV>& exactLevelInterfaceNetwork = interfaceNetworks.back()->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 j=0; j<exactLocalInterfaceAssemblers.size(); j++) {
-        //const double pen = penaltyFactor;
-        const int k = minCantorResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(j);
-        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
-        exactLocalInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
-    }
-
-    exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
-    exactGlobalAssembler.solve();
-
-    solutions[exactLevelIdx] = exactGlobalAssembler.getSol();
-    approxErrors[maxLevelIdx] = 0;
-
-    const DGBasis& exactBasis = exactGlobalAssembler.basis();
-
-    EnergyNorm<MatrixType, VectorType> energyNorm(exactGlobalAssembler.matrix());
-
-    for (int i=maxLevelIdx; i>=0; i--) {
-        const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetworks[i]->levelInterfaceNetwork(maxLevelIdx);
-        LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> globalAssembler(levelInterfaceNetwork);
-        LocalOscAssembler localAssembler(oscDataHierarchy[maxLevelIdx]->data(), oscDataHierarchy[maxLevelIdx]->mapper());
-
-        std::vector<std::shared_ptr<LocalInterfaceAssembler>> localInterfaceAssemblers(levelInterfaceNetwork.size());
-        for (size_t j=0; j<localInterfaceAssemblers.size(); j++) {
-            //const double pen = penaltyFactor;
-            const int k = minCantorResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j);
-            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
-            localInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
-        }
-
-        globalAssembler.assemble(localAssembler, localInterfaceAssemblers, l2FunctionalAssembler);
-        globalAssembler.solve();
-
-        const VectorType& sol = globalAssembler.getSol();
-        const DGBasis& basis = globalAssembler.basis();
-
-        DGMGTransfer<DGBasis> levelTransfer(basis, exactBasis);
-        levelTransfer.prolong(sol, solutions[i]);
-
-        if (i==maxLevelIdx) {
-            VectorType discErrorV(solutions[exactLevelIdx]);
-            discErrorV -= solutions[maxLevelIdx];
-            //approxErrors[i] = approxError.two_norm();
-            discError= energyNorm.operator ()(discErrorV);
-        } else {
-            VectorType approxError(solutions[maxLevelIdx]);
-            approxError -= solutions[i];
-            //approxErrors[i] = approxError.two_norm();
-            approxErrors[i] = energyNorm.operator ()(approxError);
-        }
-    }
-
-    // print results
-    std::cout << std::endl << std::endl;
-    std::cout << "cantorLevel  h     approxError"<< std::endl;
-    std::cout << "------------------------------" << std::endl;
-    for (size_t i=0; i<approxErrors.size(); i++) {
-        std::cout << (minCantorResolution+i) << "            2^-" << maxResolution << "  " << approxErrors[i] << std::endl;
-    }
-
-    std::cout << std::endl << std::endl;
-    std::cout << "exactResolution: 2^-" << exactResolution << std::endl;
-    std::cout << "discError: " << discError << std::endl;
-
-    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/CMakeLists.txt

 add_executable("geofaultnetwork" geofaultnetwork.cc)
 target_link_dune_default_libraries("geofaultnetwork")
 
-add_executable("geoconvergence" geoconvergence.cc)
-target_link_dune_default_libraries("geoconvergence")
-
 add_custom_target(geofaultnetworks_srcs SOURCES
 geofaultnetwork.parset
-geoconvergence.parset
 )
 
 

src/geofaultnetworks/geoconvergence.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/geofaultfactory.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::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("/home/mi/podlesny/software/dune/faultnetworks/src/geofaultnetworks/geoconvergence.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 minResolution = problemParameters.get<int>("minResolution");
-    const int maxResolution = problemParameters.get<int>("maxResolution");
-    const int exactResolution = problemParameters.get<int>("exactResolution");
-    const int faultCount = problemParameters.get<int>("faultCount");
-
-    const double c = problemParameters.get<double>("penaltyFactor");
-    const int patchDepth = problemParameters.get<int>("patchDepth");
-
-    const int maxIterations = problemParameters.get<int>("maxIterations");
-    const double solverTolerance = problemParameters.get<double>("tol");
-
-
-    const int maxGridN = std::pow(2, maxResolution);
-
-
-    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minResolution) + "_" + std::to_string(maxResolution) + "_" + std::to_string(exactResolution) + "_" + std::to_string(minCantorResolution) +  "_" + std::to_string((int) c) + ".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::round(std::log2(oscGridN)) - minResolution;
-
-    if (oscGridN>maxGridN)
-        DUNE_THROW(Dune::Exception, "Provided oscData too large!");
-
-    std::cout << "OscData file read successfully!" << std::endl << std::endl;
-
-    const int maxLevelIdx = maxResolution - minResolution;
-    const int exactLevelIdx = exactResolution - minResolution;
-    std::vector<int> faultToLevel(faultCount);
-
-    for (size_t i=0; i<faultToLevel.size(); ) {
-        for (int j=maxLevelIdx; j>0 && i<faultToLevel.size(); j--) {
-            faultToLevel[i] = j;
-            i++;
-        }
-
-        for (int j=0; j<maxLevelIdx && i<faultToLevel.size(); j++) {
-            faultToLevel[i] = j;
-            i++;
-        }
-    }
-
-    // build multilevel geological fault network
-    GeoFaultFactory<GridType> faultFactory(faultToLevel, minResolution, exactLevelIdx, false);
-
-    const GridType& grid = faultFactory.grid();
-    const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
-
-    for (size_t i=0; i<interfaceNetwork.size(); i++) {
-        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
-        networkWriter.write(interfaceNetwork.levelInterfaceNetwork(i), true);
-    }
-
-
-    GeoFaultFactory<GridType> faultFactory(minResolution, exactLevelIdx, minCantorResolution, maxCantorLevels);
-    const std::vector<std::shared_ptr<InterfaceNetwork<GridType>>>& interfaceNetworks = faultFactory.interfaceNetworks();
-
-    /*
-    for (size_t i=0; i<interfaceNetworks[1]->size(); i++) {
-        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
-        networkWriter.write(interfaceNetworks[1]->levelInterfaceNetwork(i), true);
-    }*/
-
-    std::cout << "Interface networks built successfully!" << std::endl << std::endl;
-
-    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 direct solutions                    ---
-    // ----------------------------------------------------------------
-    std::cout << std::endl;
-    std::cout << "Computing direct solutions!" << std::endl;
-
-    std::vector<VectorType> solutions;
-    solutions.resize(exactLevelIdx+1);
-
-    std::vector<double> approxErrors;
-    approxErrors.resize(maxLevelIdx+1);
-
-    double discError = NAN;
-
-    OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
-    L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);
-
-    const LevelInterfaceNetwork<GV>& exactLevelInterfaceNetwork = interfaceNetworks.back()->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 j=0; j<exactLocalInterfaceAssemblers.size(); j++) {
-        //const double pen = penaltyFactor;
-        const int k = minCantorResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(j);
-        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
-        exactLocalInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
-    }
-
-    exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
-    exactGlobalAssembler.solve();
-
-    solutions[exactLevelIdx] = exactGlobalAssembler.getSol();
-    approxErrors[maxLevelIdx] = 0;
-
-    const DGBasis& exactBasis = exactGlobalAssembler.basis();
-
-    EnergyNorm<MatrixType, VectorType> energyNorm(exactGlobalAssembler.matrix());
-
-    for (int i=maxLevelIdx; i>=0; i--) {
-        const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetworks[i]->levelInterfaceNetwork(maxLevelIdx);
-        LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> globalAssembler(levelInterfaceNetwork);
-        LocalOscAssembler localAssembler(oscDataHierarchy[maxLevelIdx]->data(), oscDataHierarchy[maxLevelIdx]->mapper());
-
-        std::vector<std::shared_ptr<LocalInterfaceAssembler>> localInterfaceAssemblers(levelInterfaceNetwork.size());
-        for (size_t j=0; j<localInterfaceAssemblers.size(); j++) {
-            //const double pen = penaltyFactor;
-            const int k = minCantorResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j);
-            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
-            localInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
-        }
-
-        globalAssembler.assemble(localAssembler, localInterfaceAssemblers, l2FunctionalAssembler);
-        globalAssembler.solve();
-
-        const VectorType& sol = globalAssembler.getSol();
-        const DGBasis& basis = globalAssembler.basis();
-
-        DGMGTransfer<DGBasis> levelTransfer(basis, exactBasis);
-        levelTransfer.prolong(sol, solutions[i]);
-
-        if (i==maxLevelIdx) {
-            VectorType discErrorV(solutions[exactLevelIdx]);
-            discErrorV -= solutions[maxLevelIdx];
-            //approxErrors[i] = approxError.two_norm();
-            discError= energyNorm.operator ()(discErrorV);
-        } else {
-            VectorType approxError(solutions[maxLevelIdx]);
-            approxError -= solutions[i];
-            //approxErrors[i] = approxError.two_norm();
-            approxErrors[i] = energyNorm.operator ()(approxError);
-        }
-    }
-
-    // print results
-    std::cout << std::endl << std::endl;
-    std::cout << "cantorLevel  h     approxError"<< std::endl;
-    std::cout << "------------------------------" << std::endl;
-    for (size_t i=0; i<approxErrors.size(); i++) {
-        std::cout << (minCantorResolution+i) << "            2^-" << maxResolution << "  " << approxErrors[i] << std::endl;
-    }
-
-    std::cout << std::endl << std::endl;
-    std::cout << "exactResolution: 2^-" << exactResolution << std::endl;
-    std::cout << "discError: " << discError << std::endl;
-
-    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/geoconvergence.parset

-path        = /storage/mi/podlesny/data/osccoeff/
-resultPath  = /home/mi/podlesny/results/faultnetworks/geofaultnetworks/geoconvergence/
-
-
-[problem0]
-oscDataFile         = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-minResolution = 0
-maxResolution = 9
-exactResolution = 10
-faultCount = 30
-
-penaltyFactor = 1
-patchDepth = 1
-
-# cg solver
-tol = 1e-12
-maxIterations = 8
-useExactSol = 0
-
-###########################################

src/uniformfaultnetworks/2faultscontcoarse2level.parset

-path        = /home/mi/podlesny/data/osccoeff/
-resultPath  = /home/mi/podlesny/results/faultnetworks/uniformfaultnetworks/2level/2faults/contCoarse/
-
-# 2
-############################################
-
-
-[problem0]
-oscDataFile         = oscDataLaplace64.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 9
-exactResolution = 10
-faultDistance = 2
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-
-
-[problem1]
-oscDataFile         = oscDataLaplace64.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 9
-exactResolution = 10
-faultDistance = 5
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-
-[problem2]
-oscDataFile         = oscDataLaplace64.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 9
-exactResolution = 10
-faultDistance = 15
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-
-
-[problem3]
-oscDataFile         = oscDataLaplace64.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 9
-exactResolution = 10
-faultDistance = 25
-
-# cg solver
-tol = 1e-12
-maxIterations = 20

src/uniformfaultnetworks/2faultsdiscontcoarse2level.parset

-path        = /home/mi/podlesny/data/osccoeff/
-resultPath  = /home/mi/podlesny/results/faultnetworks/uniformfaultnetworks/2level/2faults/discontCoarse/
-
-# 2
-############################################
-
-
-
-[problem0]
-oscDataFile         = oscDataLaplace64.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 4
-fineResolution = 7
-exactResolution = 10
-faultDistance = 3
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-
-
-[problem1]
-oscDataFile         = oscDataLaplace64.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 5
-fineResolution = 7
-exactResolution = 10
-faultDistance = 6
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-
-
-[problem2]
-oscDataFile         = oscDataLaplace64.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 7
-exactResolution = 10
-faultDistance = 12
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-
-
-[problem3]
-oscDataFile         = oscDataLaplace64.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 4
-fineResolution = 8
-exactResolution = 10
-faultDistance = 3
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-
-[problem4]
-oscDataFile         = oscDataLaplace64.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 4
-fineResolution = 9
-exactResolution = 10
-faultDistance = 3
-
-# cg solver
-tol = 1e-12
-maxIterations = 20

src/uniformfaultnetworks/CMakeLists.txt

-#add_executable("faultnetworks" faultnetworks.cc)
-#target_link_dune_default_libraries("faultnetworks")
-
-#add_executable("discontcoarse2level" discontcoarse2level.cc)
-#target_link_dune_default_libraries("discontcoarse2level")
-
-#add_executable("2faultscontcoarse2level" 2faultscontcoarse2level.cc)
-#target_link_dune_default_libraries("2faultscontcoarse2level")
-
-#add_executable("2faultsdiscontcoarse2level" 2faultsdiscontcoarse2level.cc)
-#target_link_dune_default_libraries("2faultsdiscontcoarse2level")
-
-#add_executable("multilevelfaultnetworks" multilevelfaultnetworks.cc)
-#target_link_dune_default_libraries("multilevelfaultnetworks")
-
-add_executable("multilevelpatches" multilevelpatches.cc)
-target_link_dune_default_libraries("multilevelpatches")
-
-
-add_custom_target(uniformfaultnetworks_srcs SOURCES
-# faultnetworks.parset
-# discontcoarse2level.parset
-# 2faultscontcoarse2level.parset
-# 2faultsdiscontcoarse2level.parset
-# multilevelfaultnetworks.parset
- multilevelpatches.parset
-)
-
-

src/uniformfaultnetworks/discontcoarse2level.parset

-path        = /home/mi/podlesny/data/osccoeff/
-resultPath  = /home/mi/podlesny/results/faultnetworks/uniformfaultnetworks/2level/discontCoarse/
-
-# 2
-############################################
-
-
-
-[problem0]
-oscDataFile         = oscDataLaplace64.mat
-
-#amplitude = 0
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 9
-exactResolution = 10
-faultCount = 63
-
-# cg solver
-tol = 1e-12
-maxIterations = 20

src/uniformfaultnetworks/faultnetworks.parset

-path        = /home/mi/podlesny/data/osccoeff/
-resultPath  = /home/mi/podlesny/results/faultnetworks/uniformfaultnetworks/2level/contCoarse/laplace/
-
-# 2
-############################################
-
-
-[problem0]
-oscDataFile         = oscDataLaplace64.mat
-exactSolFile        = oscDataLaplace64.vec
-isLaplace           = 0
-
-#amplitude = 0
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 7
-exactResolution = 10
-faultCount = 1
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-useExactSol = 0
-
-[problem1]
-oscDataFile         = oscDataLaplace64.mat
-exactSolFile        = oscDataLaplace64.vec
-isLaplace           = 0
-
-#amplitude = 0
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 7
-exactResolution = 10
-faultCount = 2
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-useExactSol = 0
-
-[problem2]
-oscDataFile         = oscDataLaplace64.mat
-exactSolFile        = oscDataLaplace64.vec
-isLaplace           = 0
-
-#amplitude = 0
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 8
-exactResolution = 10
-faultCount = 1
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-useExactSol = 0
-
-[problem3]
-oscDataFile         = oscDataLaplace64.mat
-exactSolFile        = oscDataLaplace64.vec
-isLaplace           = 0
-
-#amplitude = 0
-
-# level resolution in 2^(-...)
-coarseResolution = 6
-fineResolution = 8
-exactResolution = 10
-faultCount = 2
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-useExactSol = 0
-

src/uniformfaultnetworks/multilevelfaultnetworks.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
-*/
-
-#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/istl/matrix.hh>
-#include <dune/common/timer.hh>
-
-// dune-istl includes
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/preconditioners.hh>
-
-// dune-fufem includes
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/assemblers/assembler.hh>
-#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-
-// dune-grid includes
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/io/file/vtk/vtkwriter.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-#include <dune/grid/common/gridfactory.hh>
-#include <dune/grid/common/mcmgmapper.hh>
-
-// dune-solver includes
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/norms/twonorm.hh>
-//#include <dune/solvers/norms/onenorm.hh>
-#include <dune/faultnetworks/compressedmultigridtransfer.hh>
-
-#include <dune/faultnetworks/utils/matrixreader.hh>
-#include <dune/faultnetworks/utils/vectorreader.hh>
-#include <dune/faultnetworks/oscrhs.hh>
-
-#include <dune/faultnetworks/preconditioners/multilevelfaultpreconditioner.hh>
-
-#include <dune/faultnetworks/solvers/osccgsolver.hh>
-
-#include <dune/faultnetworks/faultp1nodalbasis.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/faultinterface.hh>
-#include <dune/faultnetworks/assemblers/interfaceoperatorassembler.hh>
-#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
-#include <dune/faultnetworks/assemblers/osclocalassembler.hh>
-#include <dune/faultnetworks/faultfactories/oscunitcube.hh>
-#include <dune/faultnetworks/faultfactories/multileveluniformfaultfactory.hh>
-#include <dune/faultnetworks/dgmgtransfer.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-#include <dune/faultnetworks/oscdata.hh>
-#include <dune/faultnetworks/oscdatahierarchy.hh>
-
-#include <dune/faultnetworks/levelinterfacenetworkproblem.hh>
-
-#include <dune/istl/superlu.hh>
-#include <dune/istl/umfpack.hh>
-
-
-//using namespace Dune;
-
-
-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;
-//    typedef OscUnitCube<GridType, 2> GridOb;
-    #else
-    #error No UG!
-    #endif
-
-
-    typedef typename GridType::LevelGridView GV;
-
-//    typedef P1NodalBasis<GV, double> Basis;
-//    typedef typename Basis::LocalFiniteElement FE;
-
-    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("/home/mi/podlesny/software/dune/faultnetworks/src/uniformfaultnetworks/multilevelfaultnetworks.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 int faultCount = problemParameters.get<int>("faultCount");
-
-    const int maxIterations = problemParameters.get<int>("maxIterations");
-    const double solverTolerance = problemParameters.get<double>("tol");
-    //const bool useExactSol = problemParameters.get<bool>("useExactSol");
-
-    //const int coarseGridN = std::pow(2, coarseResolution);
-    const int fineGridN = std::pow(2, fineResolution);
-
-
-    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(coarseResolution) + "_" + std::to_string(fineResolution) + "_" + std::to_string(exactResolution) + "_" + std::to_string(faultCount) + ".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::round(std::log2(oscGridN)) - coarseResolution;
-
-    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;
-
-    std::vector<int> faultToLevel(faultCount);
-    for (size_t i=0; i<faultToLevel.size(); ) {
-        for (int j=fineLevelIdx; j>0 && i<faultToLevel.size(); j--) {
-            faultToLevel[i] = j;
-            i++;
-        }
-
-        for (int j=0; j<fineLevelIdx && i<faultToLevel.size(); j++) {
-            faultToLevel[i] = j;
-            i++;
-        }
-    }
-
-    print(faultToLevel, "faultToLevel: ");
-
-    MultilevelUniformFaultFactory<GridType> faultFactory(faultCount, faultToLevel, coarseResolution, exactLevelIdx);
-    faultFactory.prolongToAll();
-    const GridType& grid = faultFactory.grid();
-    const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
-
-    for (size_t i=0; i<interfaceNetwork.size(); i++) {
-        std::cout << "Level: " << i << " faults: " << interfaceNetwork.levelInterfaceNetwork(i).size() << std::endl;
-    }
-
-    OscData<GridType, OscMatrixType> oscData(grid, oscGridLevelIdx);
-    oscData.set(matrix);
-    OscDataHierarchy<GridType, OscMatrixType> oscDataHierarchy(oscData);
-    
-    Timer totalTimer;
-    totalTimer.reset();
-    totalTimer.start();
-
-    // ----------------------------------------------------------------
-    // ---              compute initial iterate                     ---
-    // ----------------------------------------------------------------
-    std::cout << std::endl;
-    std::cout << "Computing initial iterate!" << std::endl;
-
-    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> coarseGlobalAssembler(interfaceNetwork.levelInterfaceNetwork(0));
-    LocalOscAssembler coarseLocalAssembler(oscDataHierarchy[0]->data(), oscDataHierarchy[0]->mapper());
-
-    OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
-    L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);
-
-    std::shared_ptr<LocalInterfaceAssembler> coarseLocalInterfaceAssemblerPtr = std::make_shared<LocalInterfaceAssembler>(B);
-    std::vector<std::shared_ptr<LocalInterfaceAssembler>> coarseLocalInterfaceAssemblers(interfaceNetwork.levelInterfaceNetwork(0).size(), coarseLocalInterfaceAssemblerPtr);
-
-    coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
-    coarseGlobalAssembler.solve();
-
-    const VectorType& coarseSol = coarseGlobalAssembler.getSol();
-    const DGBasis& coarseBasis = coarseGlobalAssembler.basis();
-
-    // ----------------------------------------------------------------
-    // ---              compute exact solution                      ---
-    // ----------------------------------------------------------------
-    std::cout << std::endl;
-    std::cout << "Computing exact solution!" << std::endl;
- 
-/*    const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(exactLevelIdx);
-
-    // build leaf view of all faults
-    LevelInterfaceNetwork<GV> leafInterfaceNetwork(levelInterfaceNetwork.grid(), levelInterfaceNetwork.level(), levelInterfaceNetwork.referencePoint());
-
-    for (int i=0; i<=exactLevelIdx; i++) {
-        interfaceNetwork.prolongLevelInterfaces(i, leafInterfaceNetwork);
-    }*/
-
-    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++) {
-        //exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, 1+3*exactLevelInterfaceNetwork.getInterface(i)->coarsestFather()->level());
-        exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B);
-    }
-
-    exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
-    exactGlobalAssembler.solve();
-
-    const VectorType& exactSol = exactGlobalAssembler.getSol();
-    const 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++) {
-        //fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, 1+3*fineLevelInterfaceNetwork.getInterface(i)->coarsestFather()->level());
-        fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B);
-    }
-
-    fineGlobalAssembler.assemble(fineLocalAssembler, fineLocalInterfaceAssemblers, l2FunctionalAssembler);
-    fineGlobalAssembler.solve();
-
-    const VectorType& fineSol = fineGlobalAssembler.getSol();
-    const VectorType& fineRhs = fineGlobalAssembler.rhs();
-    const DGBasis& fineBasis = fineGlobalAssembler.basis();
-
-    EnergyNorm<MatrixType, VectorType> fineEnergyNorm(fineGlobalAssembler.matrix());
-
-    // ---------------------------------------------------------------
-    // multilevel transfer
-    // ---------------------------------------------------------------
-    DGMGTransfer<DGBasis> coarseFineTransfer(coarseBasis, fineBasis);
-
-
-
-    // TODO: remove after debugging
-    /*typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> P1FECache;
-    typedef typename P1FECache::FiniteElementType FEType;
-    P1FECache p1FECache;
-    typedef typename GV::template Codim<0>::Iterator ElementIterator;
-
-    ElementIterator fIt    = fineGridView.template begin<0>();
-    ElementIterator fEndIt = fineGridView.template end<0>();
-
-    for (; fIt != fEndIt; ++fIt) {
-        int elemIdx = fineMapper.index(*fIt);
-
-        const Dune::ReferenceElement<double,dim>& fineRefElement
-            = Dune::ReferenceElements<double, dim>::general(fIt->type());
-
-        // Get local finite element
-        const FEType& fineBaseSet = p1FECache.get(fIt->type());
-
-        const size_t numFineBaseFct = fineBaseSet.localBasis().size();
-
-        std::vector<Dune::FieldVector<double,dim> > local(numFineBaseFct);
-
-        for (size_t i=0; i<numFineBaseFct; i++)
-        {
-            const Dune::LocalKey& iLocalKey = fineBaseSet.localCoefficients().localKey(i);
-            local[i] = fineRefElement.position(iLocalKey.subEntity(), iLocalKey.codim());
-            int globalFine = fineBasis.index(*fIt, iLocalKey.subEntity());
-
-            std::cout << "ElementIdx: " << elemIdx << "   " << globalFine << std::endl;
-        }
-    } */
-    // end remove
-
-
-    std::cout << "Setting up the preconditioner!" << std::endl;
-   /* BitVector activeLevels(fineLevelIdx+1, false);
-    activeLevels[fineLevelIdx][0] = true;
- */
-
-    BitVector activeLevels(fineLevelIdx+1, true);
-    //activeLevels[0][0] = false;
-    //activeLevels[1][0] = false;
-    //activeLevels[2][0] = false;
-    //activeLevels[3][0] = false;*/
-
-    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++) {
-            //levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, 1+3*levelInterfaceNetwork.getInterface(j)->coarsestFather()->level());
-            levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B);
-        }
-    }
-
-    typedef MultilevelFaultPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelFaultPreconditioner;
-    typedef typename MultilevelFaultPreconditioner::LevelFaultPreconditionerType LevelFaultPreconditioner;
-    MultilevelFaultPreconditioner::Mode preconditionerMode = MultilevelFaultPreconditioner::Mode::additive;
-
-    MultilevelFaultPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
-    for (size_t i=0; i<preconditioner.size(); i++) {
-        LevelFaultPreconditioner& levelFaultPreconditioner = *preconditioner.levelFaultPreconditioner(i);
-        Dune::BitSetVector<1> activeFaults(levelFaultPreconditioner.levelInterfaceNetwork().size(), true);
-
-        levelFaultPreconditioner.setActiveFaults(activeFaults);
-        levelFaultPreconditioner.setPatchDepth(0);
-        levelFaultPreconditioner.setBoundaryMode(LevelFaultPreconditioner::BoundaryMode::homogeneous);
-
-        levelFaultPreconditioner.build();
-    }
-
-    std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
-    std::cout << std::endl << std::endl;
-    
-    /*Dune::BlockVector<Dune::FieldVector<double,1>> initialError(fineX.size());
-    initialError = initialX;
-    initialError -= fineX;
-    std::cout << "Initial error in energy norm: " << energyNorm.operator()(initialError) << std::endl;
-    std::cout << "Initial error in two norm: " << initialError.two_norm() << std::endl << std::endl;*/
-
-
-    VectorType initialX;
-    coarseFineTransfer.prolong(coarseSol, initialX);
-
-    VectorType initialXCopy(initialX);
-    VectorType fineRhsCopy(fineRhs);
-    VectorType rhs(fineRhs);
-
-    //print(initialX, "initialX");
-    //print(fineRhs, "fineRhs");
-    //print(exactSol, "exactSol");
-    //print(fineX, "fineX");
-    
-    /* Basis fineContBasis(fineGridView);
-    std::cout << "fineContBasis.size(): " << fineContBasis.size() << std::endl;
-    DGMGTransfer<Basis> contMGTransfer(coarseBasis, fineContBasis);
-
-    Basis exactContBasis(exactGridView);
-    std::cout << "exactContBasis.size(): " << exactContBasis.size() << std::endl; */
-
-    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)
-    solver.check();
-    solver.preprocess();
-    solver.solve();
-
-    VectorType finalError(fineSol.size());
-    finalError = initialX;
-    finalError -= fineSol;
-    std::cout << "Final error in energy norm: " << fineEnergyNorm.operator()(finalError) << std::endl;
-    std::cout << "Final error in two norm: " << finalError.two_norm() << std::endl << std::endl;
-
-    std::cout << std::endl;
-    std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl; 
-
-    // regular unpreconditioned cg solver for comparisson
-/*
-    std::cout << "Setup complete, starting cg iteration!" << std::endl;
-    std::cout << std::endl;
-    
-    // solve
-    OscCGSolver<MatrixType, VectorType, DGBasis >
-    regSolver(fineBasis, &fineGlobalAssembler.matrix(), &initialXCopy, &fineRhsCopy, nullptr,
-                    maxIterations, solverTolerance, &exactEnergyNorm,
-                    Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true);
-    regSolver.check();
-    regSolver.preprocess();
-    regSolver.solve();
-*/
-    // TODO: use once debugging finished
-    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/uniformfaultnetworks/multilevelfaultnetworks.parset

-path        = /storage/mi/podlesny/data/osccoeff/
-resultPath  = /home/mi/podlesny/results/faultnetworks/uniformfaultnetworks/finelevel/
-
-
-[problem0]
-oscDataFile         = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 4
-fineResolution = 7
-exactResolution = 8
-faultCount = 15
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-useExactSol = 0

src/uniformfaultnetworks/multilevelpatches.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
-*/
-
-#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/bcrsmatrix.hh>
-#include <dune/istl/matrix.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/matrixreader.hh>
-#include <dune/faultnetworks/utils/vectorreader.hh>
-#include <dune/faultnetworks/utils/levelinterfacenetworkwriter.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-#include <dune/faultnetworks/assemblers/interfaceoperatorassembler.hh>
-#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
-#include <dune/faultnetworks/assemblers/osclocalassembler.hh>
-
-#include <dune/faultnetworks/solvers/osccgsolver.hh>
-
-#include <dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh>
-
-#include <dune/faultnetworks/faultp1nodalbasis.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/faultinterface.hh>
-#include <dune/faultnetworks/levelinterfacenetworkproblem.hh>
-
-#include <dune/faultnetworks/dgmgtransfer.hh>
-#include <dune/faultnetworks/compressedmultigridtransfer.hh>
-
-#include <dune/faultnetworks/oscdata.hh>
-#include <dune/faultnetworks/oscdatahierarchy.hh>
-#include <dune/faultnetworks/oscrhs.hh>
-
-#include <dune/faultnetworks/faultfactories/multileveluniformfaultfactory.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("/home/mi/podlesny/software/dune/faultnetworks/src/uniformfaultnetworks/multilevelpatches.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 int faultCount = problemParameters.get<int>("faultCount");
-
-    const double penaltyFactor = problemParameters.get<double>("penaltyFactor");
-    const int patchDepth = problemParameters.get<int>("patchDepth");
-
-    const int maxIterations = problemParameters.get<int>("maxIterations");
-    const double solverTolerance = problemParameters.get<double>("tol");
-    //const bool useExactSol = problemParameters.get<bool>("useExactSol");
-
-    //const int coarseGridN = std::pow(2, coarseResolution);
-    const int fineGridN = std::pow(2, fineResolution);
-
-
-    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(coarseResolution) + "_" + std::to_string(fineResolution) + "_" + std::to_string(exactResolution) + "_" + std::to_string(faultCount) + ".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::round(std::log2(oscGridN)) - coarseResolution;
-
-    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;
-
-    std::vector<int> faultToLevel(faultCount);
-    for (size_t i=0; i<faultToLevel.size(); ) {
-        for (int j=fineLevelIdx; j>0 && i<faultToLevel.size(); j--) {
-            faultToLevel[i] = j;
-            i++;
-        }
-
-        for (int j=0; j<fineLevelIdx && i<faultToLevel.size(); j++) {
-            faultToLevel[i] = j;
-            i++;
-        }
-    }
-
-    print(faultToLevel, "faultToLevel: ");
-
-    MultilevelUniformFaultFactory<GridType> faultFactory(faultCount, faultToLevel, coarseResolution, exactLevelIdx);
-    faultFactory.prolongToAll();
-
-    const GridType& grid = faultFactory.grid();
-    const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
-
-    int writeLevel = fineLevelIdx;
-    LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
-    networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel));
-
-    OscData<GridType, OscMatrixType> oscData(grid, oscGridLevelIdx);
-    oscData.set(matrix);
-    OscDataHierarchy<GridType, OscMatrixType> oscDataHierarchy(oscData);
-    
-    Timer totalTimer;
-    totalTimer.reset();
-    totalTimer.start();
-
-
-    // ----------------------------------------------------------------
-    // ---              compute initial iterate                     ---
-    // ----------------------------------------------------------------
-    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());
-
-    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 double pen = penaltyFactor;
-        const double pen = penaltyFactor*std::pow(2.0, coarseResolution);
-        coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
-    }
-
-    coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
-    coarseGlobalAssembler.solve();
-
-    const VectorType& coarseSol = coarseGlobalAssembler.getSol();
-    const DGBasis& coarseBasis = coarseGlobalAssembler.basis();
-
-    // ----------------------------------------------------------------
-    // ---              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 double pen = penaltyFactor;
-        const double pen = penaltyFactor*std::pow(2.0, coarseResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(i));
-        exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
-    }
-
-    exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
-    exactGlobalAssembler.solve();
-
-    const VectorType& exactSol = exactGlobalAssembler.getSol();
-    const 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 double pen = penaltyFactor;
-        const double pen = penaltyFactor*std::pow(2.0, coarseResolution + fineLevelInterfaceNetwork.getIntersectionsLevels().at(i));
-        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();
-    const 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);
-    //activeLevels[0][0] = true;
-    //activeLevels[1][0] = false;
-    //activeLevels[2][0] = false;
-    //activeLevels[fineLevelIdx][0] = 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 double pen = penaltyFactor;
-            const double pen = penaltyFactor*std::pow(2.0, coarseResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j));
-            levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
-        }
-    }
-
-    typedef MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelPatchPreconditioner;
-    typedef typename MultilevelPatchPreconditioner::LevelPatchPreconditionerType LevelPatchPreconditioner;
-    //typedef typename MultilevelPatchPreconditioner::LevelGlobalPreconditionerType LevelGlobalPreconditioner;
-    MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::additive;
-
-    MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
-    for (size_t i=0; i<preconditioner.size()-1; i++) {
-        LevelPatchPreconditioner& 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;
-    
-    /*Dune::BlockVector<Dune::FieldVector<double,1>> initialError(fineX.size());
-    initialError = initialX;
-    initialError -= fineX;
-    std::cout << "Initial error in energy norm: " << energyNorm.operator()(initialError) << std::endl;
-    std::cout << "Initial error in two norm: " << initialError.two_norm() << std::endl << std::endl;*/
-
-
-    VectorType initialX;
-    coarseFineTransfer.prolong(coarseSol, initialX);
-
-    VectorType initialXCopy(initialX);
-    VectorType fineRhsCopy(fineRhs);
-    VectorType rhs(fineRhs);
-
-    //print(initialX, "initialX");
-    //print(fineRhs, "fineRhs");
-    //print(exactSol, "exactSol");
-    //print(fineSol, "fineSol");
-    
-    /* Basis fineContBasis(fineGridView);
-    std::cout << "fineContBasis.size(): " << fineContBasis.size() << std::endl;
-    DGMGTransfer<Basis> contMGTransfer(coarseBasis, fineContBasis);
-
-    Basis exactContBasis(exactGridView);
-    std::cout << "exactContBasis.size(): " << exactContBasis.size() << std::endl; */
-
-    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)
-    solver.check();
-    solver.preprocess();
-    solver.solve();
-
-
-    /*
-    VectorType finalError(fineSol.size());
-    finalError = initialX;
-    finalError -= fineSol;
-    std::cout << "Final error in energy norm: " << fineEnergyNorm.operator()(finalError) << std::endl;
-    std::cout << "Final error in two norm: " << finalError.two_norm() << std::endl << std::endl;*/
-
-    std::cout << std::endl;
-    std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl; 
-
-    // regular unpreconditioned cg solver for comparisson
-/*
-    std::cout << "Setup complete, starting cg iteration!" << std::endl;
-    std::cout << std::endl;
-    
-    // solve
-    OscCGSolver<MatrixType, VectorType, DGBasis >
-    regSolver(fineBasis, &fineGlobalAssembler.matrix(), &initialXCopy, &fineRhsCopy, nullptr,
-                    maxIterations, solverTolerance, &exactEnergyNorm,
-                    Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true);
-    regSolver.check();
-    regSolver.preprocess();
-    regSolver.solve();
-*/
-    // TODO: use once debugging finished
-    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/uniformfaultnetworks/multilevelpatches.parset

-path        = /storage/mi/podlesny/data/osccoeff/
-resultPath  = /home/mi/podlesny/results/faultnetworks/uniformfaultnetworks/multilevelpatches/testing/
-
-
-[problem0]
-oscDataFile         = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 4
-fineResolution = 5
-exactResolution = 6
-faultCount = 1
-
-penaltyFactor = 10
-patchDepth = 1
-
-# cg solver
-tol = 1e-12
-maxIterations = 20
-useExactSol = 0