.

dune/faultnetworks/faultfactories/CMakeLists.txt

@@ -4,6 +4,7 @@ add_custom_target(faultnetworks_faultfactories_src SOURCES
   geofaultfactory.hh
   cantorconvergencefaultfactory.hh
   cantorfaultfactory.hh
+  rockfaultfactory.hh
   sparsecantorfaultfactory.hh
 )
 
@@ -12,5 +13,6 @@ install(FILES
   geofaultfactory.hh
   cantorconvergencefaultfactory.hh
   cantorfaultfactory.hh
+  rockfaultfactory.hh
   sparsecantorfaultfactory.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/faultfactories)

dune/faultnetworks/faultfactories/rockfaultfactory.hh

+#ifndef ROCK_FAULT_FACTORY_HH
+#define ROCK_FAULT_FACTORY_HH
+
+#include <math.h>
+#include <queue>
+
+#include <dune/faultnetworks/facehierarchy.hh>
+#include <dune/faultnetworks/utils/debugutils.hh>
+#include <dune/faultnetworks/faultfactories/oscunitcube.hh>
+#include <dune/faultnetworks/levelinterfacenetwork.hh>
+#include <dune/faultnetworks/interfacenetwork.hh>
+#include <dune/faultnetworks/faultinterface.hh>
+#include <dune/faultnetworks/hierarchicleveliterator.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/boundaryiterator.hh>
+
+template<class ctype = double> class Rock {
+protected:
+    int level_;
+
+    // vertex IDs
+    int left_;
+    int right_;
+    int top_;
+    int bottom_;
+
+    Dune::FieldVector<ctype, 2> center_;
+public:
+    void set(int level, int left, int right, int top, int bottom, const Dune::FieldVector<ctype, 2>& center) {
+        level_ = level;
+
+        left_ = left;
+        right_ = right;
+        top_ = top;
+        bottom_ = bottom;
+
+        center_ = center;
+    }
+
+    auto level() const {
+        return level_;
+    }
+
+    auto left() const {
+        return left_;
+    }
+
+    auto right() const {
+        return right_;
+    }
+
+    auto top() const {
+        return top_;
+    }
+
+    auto bottom() const {
+        return bottom_;
+    }
+
+    const auto& center() {
+        return center_;
+    }
+};
+
+
+template <class GridType>
+class LevelRockFaultFactory
+{
+    //! Parameter for mapper class
+    template<int dim>
+    struct FaceMapperLayout
+    {
+        bool contains (Dune::GeometryType gt)
+        {
+            return gt.dim() == dim-1;
+        }
+    };
+
+protected:
+    static const int dimworld = GridType::dimensionworld;
+    static const int dim = GridType::dimension;
+    using ctype = typename GridType::ctype;
+
+    using Coords = typename Dune::FieldVector<ctype, dimworld>;
+    using GV = typename GridType::LevelGridView;
+    using Intersection = typename GridType::LevelIntersection;
+
+    using Element = typename GridType::template Codim<0>::Entity;
+    static const int dimElement = Element::dimension;
+
+    using FaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout >;
+
+
+    using Rock = Rock<ctype>;
+
+    const std::shared_ptr<GridType> grid_;
+    const int level_;
+    const ctype resolution_;
+    const GV gridView_;
+    const double maxAngle_;
+
+    const typename GV::IndexSet& indexSet_;
+
+    FaceMapper faceMapper_;
+
+    std::vector<Intersection> faces_;
+    std::vector<Coords> vertexPositions_;
+    std::vector<std::vector<size_t>> vertexToFaces_;
+    std::vector<int> coarseToLevelVertex_;
+
+    const LevelRockFaultFactory& coarseLevelFactory_;
+    std::vector<Rock> rocks_;
+    std::vector<std::shared_ptr<FaultInterface<GV>>>& faults_;
+
+private:
+    template <typename T> int sgn(T val) {
+        return (T(0) < val) - (val < T(0));
+    }
+
+        bool intersectionAllowed(const size_t faceID, const size_t vertexID,
+                                 const std::array<size_t, 2>& centerIDs,
+                                 const std::set<int>& separatingIDs, const std::set<size_t>& faultDofs,
+                                 const Coords& direction, size_t dim) const {
+
+            const auto& intersection = faces_[faceID];
+
+            //check if "back" edge, cannot deviate from desiredOrientation more than maxAngle degrees (radian)
+            auto orientation = intersection.geometry().center();
+            orientation -= vertexPositions_[vertexID];
+            orientation /= orientation.two_norm();
+
+            if (std::acos(direction*orientation) > maxAngle_)
+                return false;
+
+            // check if intersection has separating dofs or other fault dofs
+            std::set<size_t> intersectionDofs;
+            computeIntersectionDofs(indexSet_, intersection, intersectionDofs);
+            intersectionDofs.erase(vertexID);
+
+            const size_t otherDim = (dim + 1) % 2;
+            for (const auto& isDof : intersectionDofs) {
+                const auto& vertex = vertexPositions_[isDof];
+                std::array<size_t, 2> IDs = {computeID(vertex, 0), computeID(vertex, 1)};
+
+                bool centerPassed = (sgn(direction[dim])>0) ? (centerIDs[dim]<IDs[dim]) : (centerIDs[dim]>IDs[dim]);
+                if (faultDofs.count(isDof) or separatingIDs.count(IDs[otherDim]) or centerPassed) {
+                    return false;
+                }
+            }
+
+            return true;
+        }
+
+        /*
+        ctype distance(const Intersection& isec, const Dune::FieldVector<ctype, dimworld> & vertex) const {
+            Dune::FieldVector<ctype, dimworld> vec(vertex);
+            vec += isec.center();
+
+            return isec.unitOuterNormal()*vec;
+        }*/
+
+        //works only in 2D, vertex(1) = vertex(0) - yid, intersection of line (given by vertex and derivative 1) with y axis
+        int computeID(const Coords& vertex, int dim) const {
+            return (int) (vertex[dim]*resolution_);
+        }
+
+        void computeIntersectionDofs(const Intersection& intersection, std::set<size_t>& intersectionDofs) const {
+            intersectionDofs.clear();
+
+            // loop over all vertices of the intersection
+            const auto& insideElement = intersection.inside();
+            const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
+
+            for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
+                size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
+                size_t globalIdx = indexSet_.subIndex(insideElement, idxInElement, dimElement);
+
+                intersectionDofs.insert(globalIdx);
+            }
+        }
+
+        bool isTargetBoundaryVertex(size_t vertexIdx, const std::array<size_t, 2>& targetIDs) const {
+            const auto& vertex = vertexPositions_[vertexIdx];
+            const std::array<size_t, 2> vertexIDs = {computeID(vertex, 0), computeID(vertex, 1)};
+            return (vertexIDs[0] == targetIDs[0]) and (vertexIDs[1] == targetIDs[1]);
+        }
+
+        typedef std::vector<double>::const_iterator ConstVectorIterator;
+
+        struct ordering {
+            bool operator ()(std::pair<size_t, ConstVectorIterator> const& a, std::pair<size_t, ConstVectorIterator> const& b) {
+                return *(a.second) < *(b.second);
+            }
+        };
+
+        void distanceSort(std::vector<size_t>& admissibleFaces, std::vector<double>& distances) const {
+            std::vector<std::pair<size_t, ConstVectorIterator>> order(distances.size());
+
+            size_t j = 0;
+            ConstVectorIterator it = distances.begin();
+            ConstVectorIterator itEnd = distances.end();
+            for (; it != itEnd; ++it, ++j)
+                order[j] = std::make_pair(j, it);
+
+            sort(order.begin(), order.end(), ordering());
+
+            std::vector<size_t> initialAdmissibleFaces(admissibleFaces);
+            for (size_t i=0; i<admissibleFaces.size(); i++)
+                admissibleFaces[i] = initialAdmissibleFaces[order[i].first];
+        }
+
+       /* void generateFaultSeeds(std::vector<FaultCorridor>& faultCorridors, std::vector<std::vector<size_t>>& faultSeeds) {
+
+            //std::cout << "LevelGeoFaultFactory::generateFaultSeeds() " << std::endl;
+            //std::cout << "------------------------------------- " << std::endl << std::endl;
+
+            faultSeeds.resize(faultCorridors.size());
+
+            std::map<int, std::pair<size_t,size_t>> yIDtoFaultCorridor;
+            for (size_t i=0; i<faultCorridors.size(); i++){
+                FaultCorridor& faultCorridor = faultCorridors[i];
+                const std::vector<int> & faultYs = faultCorridor.faultYs();
+
+                faultSeeds[i].resize(faultYs.size());
+                for (size_t j=0; j<faultYs.size(); j++) {
+                    yIDtoFaultCorridor[faultYs[j]] = std::make_pair(i, j);
+                }
+            }
+
+            BoundaryIterator<GV> bIt(gridView_, BoundaryIterator<GV>::begin);
+            BoundaryIterator<GV> bEnd(gridView_, BoundaryIterator<GV>::end);
+            for(; bIt!=bEnd; ++bIt) {
+
+                const Element& insideElement = bIt->inside();
+                const auto& geometry = insideElement.geometry();
+
+                // neglect "upper" part of boundary, where (y==1 and x>0) or (y>0 and x==1)
+                const auto& center = bIt->geometry().center();
+                if (center[0] + center[1] >1)
+                        continue;
+
+                const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
+
+                for (int i=0; i<refElement.size(bIt->indexInInside(), 1, dimElement); i++) {
+                    size_t idxInElement = refElement.subEntity(bIt->indexInInside(), 1, i, dimElement);
+                    size_t globalIdx = indexSet_.subIndex(insideElement, idxInElement, dimElement);
+
+                    const auto& vertex = geometry.corner(idxInElement);
+
+                    const int yID = computeYID(vertex);
+
+                    if (yIDtoFaultCorridor.count(yID)) {
+                        const std::pair<size_t, size_t>& indices = yIDtoFaultCorridor[yID];
+                        faultSeeds[indices.first][indices.second] = globalIdx;
+                    }
+                }
+            }
+
+            //std::cout << "------------------------------------- " << std::endl << std::endl;
+        }*/
+
+
+        void generateFault(std::shared_ptr<FaultInterface<GV>> fault, const size_t faultSeedIdx, const Coords& center,
+                           const size_t corridor,
+                           const std::set<size_t>& constSeparatingDofs) {
+            //std::cout << "LevelGeoFaultFactory::generateFault() " << std::endl;
+            //std::cout << "------------------------------------- " << std::endl << std::endl;
+
+            bool success = false;
+
+            std::set<size_t> separatingDofs = constSeparatingDofs;
+
+            auto direction = center;
+            direction -= vertexPositions_[faultSeedIdx];
+            direction /= direction.two_norm();
+
+            size_t dim = 0;
+            if (std::abs(direction[0]) < std::abs(direction[1])) {
+                dim = 1;
+            }
+
+            const std::array<size_t, 2> centerIDs = {computeID(center, 0), computeID(center, 1)};
+
+            std::set<size_t> separatingIDs;
+
+            auto faultSeedID = computeID(vertexPositions_[faultSeedIdx], dim);
+            separatingIDs.insert(faultSeedID + corridor);
+            separatingIDs.insert(faultSeedID - corridor);
+
+            std::vector<size_t> faultDofs(0);
+            faultDofs.push_back(faultSeedIdx);
+
+            std::vector<size_t> faultFaces;
+
+            std::map<size_t, std::vector<size_t>> vertexToAdmissibleFaces;
+
+            std::queue<size_t> vertexQueue;
+            vertexQueue.push(faultSeedIdx);
+
+            while (!vertexQueue.empty()) {
+                const size_t vertexID = vertexQueue.front();
+                vertexQueue.pop();
+
+                if (isTargetBoundaryVertex(vertexID, centerIDs)) {
+                    success = true;
+                    break;
+                }
+
+                if (vertexToAdmissibleFaces.count(vertexID)==0) {
+                    const std::vector<size_t>& faces = vertexToFaces_[vertexID];
+                    std::vector<size_t> admissibleFaces;
+
+                    for (size_t i=0; i<faces.size(); i++) {
+                        if (intersectionAllowed(faces[i], vertexID, centerIDs, separatingDofs, separatingIDs, faultDofs, direction, dim)) {
+                            admissibleFaces.push_back(faces[i]);
+                        }
+                    }
+
+                    std::vector<double> distances(admissibleFaces.size(), 0);
+                    for (size_t i=0; i<admissibleFaces.size(); i++) {
+                        Coords vec(center);
+                        vec -= faces_[admissibleFaces[i]].geometry().center();
+                        distances[i] = vec.two_norm();
+                    }
+                    distanceSort(admissibleFaces, distances);
+
+                    vertexToAdmissibleFaces[vertexID] = admissibleFaces;
+                }
+
+                std::vector<size_t>& suitableFaces = vertexToAdmissibleFaces[vertexID];
+
+                if (suitableFaces.size()==0) {
+                    faultDofs.pop_back();
+                    faultFaces.pop_back();
+                    vertexToAdmissibleFaces.erase(vertexID);
+
+                    separatingDofs.insert(vertexID);
+                    vertexQueue.push(faultDofs.back());
+                } else {
+
+                    // generate random number from (0,1)
+                    double randomNumber = ((double) std::rand() / (RAND_MAX));
+
+                    size_t randSelector = (size_t) std::abs(std::log2(randomNumber));
+                    if (randSelector >= suitableFaces.size()) {
+                        randSelector = suitableFaces.size() -1;
+                    }
+
+                    size_t nextFaceID = suitableFaces[randSelector];
+
+                    suitableFaces.erase(suitableFaces.begin()+randSelector);
+
+                    std::set<size_t> nextFaceDofs;
+                    computeIntersectionDofs(indexSet_, faces_[nextFaceID], nextFaceDofs);
+                    nextFaceDofs.erase(vertexID);
+
+                    faultDofs.push_back(*nextFaceDofs.begin());
+                    faultFaces.push_back(nextFaceID);
+
+                    size_t nextVertexID = *(nextFaceDofs.begin());
+                    vertexQueue.push(nextVertexID);
+                }
+            }
+
+            if (!success) {
+                std::cout << "Generating a fault failed: Unable to reach target boundary! This should not happend!" << std::endl;
+                DUNE_THROW(Dune::Exception, "Generating a fault failed: Unable to reach target boundary!");
+            }
+
+            for (size_t i=0; i<faultFaces.size(); i++) {
+                fault->addFace(faces_[faultFaces[i]]);
+            }
+
+            //std::cout << "------------------------------------- " << std::endl << std::endl;
+        }
+
+    void prolong(const Rock& rock, Rock& newRock) {
+        auto newLeft = coarseToLevelVertex_[rock.left()];
+        auto newRight = coarseToLevelVertex_[rock.right()];
+        auto newTop = coarseToLevelVertex_[rock.top()];
+        auto newBottom = coarseToLevelVertex_[rock.bottom()];
+
+        newRock.set(rock.level(), newLeft, newRight, newTop, newBottom, rock.center());
+    }
+
+    bool randomSplit(const Rock& rock) {
+        const auto& center = rock.center();
+
+        double res = 0.0;
+        if (center[0] < center[1]) {
+            res = 100.0*(1.0 - center[0]);
+        } else {
+            res = 100.0*(1.0 - center[1]);
+        }
+
+        double prob = (100.0 * std::rand() / (RAND_MAX + 1.0)) + 1;
+        return  res > prob;
+    }
+
+    void split(const Rock& rock, const std::set<size_t>& separatingDofs) {
+        Rock newRock;
+        prolong(rock, newRock);
+
+        bool toBeSplit = (rock.level() == level_-1) and randomSplit(rock);
+        if (!toBeSplit) {
+            rocks_.push_back(newRock);
+        } else {
+            // split rock into 4 subparts by two new intersecting faults
+            std::shared_ptr<FaultInterface<GV>> xFault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
+            generateFault(xFault, newRock, separatingDofs, 0);
+            faults_.push_back(xFault);
+
+            std::shared_ptr<FaultInterface<GV>> yFault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
+            generateFault(yFault, newRock, separatingDofs, 1);
+            faults_.push_back(yFault);
+        }
+    }
+
+    public:
+        LevelRockFaultFactory(const std::shared_ptr<GridType> grid, const int level, const ctype resolution) :
+            grid_(grid),
+            level_(level),
+            resolution_(resolution),
+            gridView_(grid_->levelGridView(level_)),
+            indexSet_(gridView_.indexSet()),
+            faceMapper_(gridView_)
+            {
+                // init faces_, vertexPositions_, vertexToFaces_
+                faces_.resize(faceMapper_.size());
+                vertexPositions_.resize(gridView_.size(dim));
+                vertexToFaces_.resize(gridView_.size(dim));
+                coarseToLevelVertex_.resize(coarseLevelFactory_.gridView().size(dim));
+
+                std::vector<bool> faceHandled(faceMapper_.size(), false);
+
+                for (const auto& elem:elements(gridView_)) {
+                    for (const auto& isect:intersections(gridView_, elem)) {
+                        //const size_t faceID = faceMapper_.index(isect);
+
+                        const int faceID = faceMapper_.subIndex(elem, isect.indexInInside(), 1);
+
+                        if (isect.boundary())
+                            continue;
+
+                        if (faceHandled[faceID])
+                            continue;
+
+                        faceHandled[faceID] = true;
+                        faces_[faceID] = isect;
+
+                        const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(elem.type());
+
+                        for (int i=0; i<refElement.size(isect.indexInInside(), 1, dimElement); i++) {
+                            //size_t vertexID = idxSet.subIndex(elem.subEntity<1>(isect), i, dimElement);
+
+                            size_t idxInElement = refElement.subEntity(isect.indexInInside(), 1, i, dimElement);
+                            const size_t vertexID = indexSet_.subIndex(elem, idxInElement, dimElement);
+
+                            vertexPositions_[vertexID] = elem.geometry().corner(idxInElement);
+                            vertexToFaces_[vertexID].push_back(faceID);
+                        }
+                    }
+                }
+
+                for (const auto& vertex : vertices(gridView_)) {
+                    size_t coarseVertexID;
+                    bool isCoarseVertex = coarseLevelFactory_.indexSet().index(vertex, coarseVertexID);
+
+                    if (isCoarseVertex) {
+                        coarseToLevelVertex_[coarseVertexID] = indexSet_.index(vertex);
+                    }
+                }
+            }
+
+        void build(const std::set<size_t>& separatingDofs) {
+            faults_.resize(0);
+            const auto& coarseRocks = coarseLevelFactory_.rocks();
+
+            for (size_t i=0; i<coarseRocks.size(); i++) {
+                split(coarseRocks[i], separatingDofs);
+            }
+        }
+
+        const auto& indexSet() const {
+            return indexSet_;
+        }
+
+        const auto& rocks() const {
+            return rocks_;
+        }
+};
+
+template <class GridType>
+class RockFaultFactory {
+
+private:
+    const int coarseResolution_;
+    const size_t maxLevel_;
+    const int coarseGridN_;
+
+    std::vector<double> levelResolutions_;
+
+    std::shared_ptr<GridType> grid_;
+    std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_;
+
+    std::vector<std::shared_ptr<LevelRockFaultFactory<GridType>>> levelRockFaultFactories_;
+
+/*typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type
+            almost_equal(ctype x, ctype y, int ulp) const
+        {
+            return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp
+                   || std::abs(x-y) < std::numeric_limits<ctype>::min();
+        }
+
+
+
+bool isSeparatingIntersection(const Intersection& intersection, const std::set<int>& separatingYIDs) const {
+    Dune::FieldVector<ctype, dimworld> faceCenter = intersection.geometry().center();
+
+    ctype yID = faceCenter[1]*coarseGridN_;
+    ctype intpart;
+
+    if (almost_equal(std::modf(yID, &intpart), 0.0, 2)) {
+        std::set<int>::iterator it = separatingYIDs.find((int) intpart);
+        if (it!=separatingYIDs.end())
+            return true;
+        else
+            return false;
+    } else
+        return false;
+}*/
+
+
+public:
+        //setup 
+    RockFaultFactory(const int coarseResolution, const size_t maxLevel, const double maxAngle = 2) :
+        coarseResolution_(coarseResolution),
+        maxLevel_(maxLevel),
+        coarseGridN_(std::pow(2, coarseResolution_)),
+        interfaceNetwork_(nullptr)
+        {   
+            using GridOb = OscUnitCube<GridType, 2>;
+            using GV = typename GridType::LevelGridView;
+
+            Dune::UGGrid<GridType::dimension>::setDefaultHeapSize(4000);
+            GridOb unitCube(coarseGridN_);
+            grid_ = unitCube.grid();
+            grid_->globalRefine(maxLevel_);
+
+            levelResolutions_.resize(maxLevel_+1);
+            levelRockFaultFactories_.resize(maxLevel_+1);
+
+            // init interface network
+            interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);
+
+            for (size_t i=0; i<=maxLevel_; i++) {
+                std::vector<std::shared_ptr<FaultInterface<GV>>> faults;
+
+                levelResolutions_[i] = std::pow(2, coarseResolution_+i);
+
+                //generate faults on level
+                levelRockFaultFactories_[i] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, i, levelResolutions_[i]);
+                levelRockFaultFactories_[i]->build(faults, maxAngle);
+
+                for (size_t j=0; j<faults.size(); j++) {
+                    interfaceNetwork_->addInterface(faults[j]);
+                }
+
+                if (i==maxLevel_)
+                    continue;
+
+                interfaceNetwork_->prolongLevel(i, i+1);
+            }
+        }
+
+    /*
+    void prolongToAll() {
+        // prolong all faults to all subsequent levels
+        for (int i=maxLevel_-1; i>=0; i--) {
+            if (interfaceNetwork_->size(i)>0) {
+                std::set<int> toLevels;
+                for (size_t j=i+1; j<=maxLevel_; j++) {
+                    toLevels.insert(j);
+                }
+
+                interfaceNetwork_->prolongLevelInterfaces(i, toLevels);
+            }
+        }
+        interfaceNetwork_->build();
+    }*/
+
+    /*void prolongToAll() {
+        // prolong all faults to all subsequent levels
+        for (int i=interfaceNetwork_->size()-1; i>=0; i--) {
+            interfaceNetwork_->prolongLevelInterfaces(i, maxLevel_);
+        }
+        interfaceNetwork_->build();
+    }*/
+
+    const GridType& grid() const {
+        return *grid_;
+	}
+	
+    /*const InterfaceNetwork<GridType>& interfaceNetwork() {
+        return *interfaceNetwork_;
+    }*/
+
+    InterfaceNetwork<GridType>& interfaceNetwork() {
+        return *interfaceNetwork_;
+    }
+};
+#endif

dune/faultnetworks/localproblem.hh

@@ -147,6 +147,10 @@ public:
 
     void updateLocalRhs(const RangeType& rhs){
         localRhs = rhs;
+        /*for (size_t i=0; i<localRhs.size(); i++) {
+            if (!boundaryDofs_[i][0])
+                localRhs[i] = rhs[i];
+        }*/
     }
 
     void updateRhsAndBoundary(const RangeType& rhs, const RangeType& boundary) {
@@ -160,6 +164,8 @@ public:
     }
 
     void solve(DomainType& x){
+        //print(boundaryDofs_, "boundaryDofs:");
+
         #if HAVE_UMFPACK
             RangeType localRhsCopy(localRhs);
             Dune::InverseOperatorResult res;

dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh

@@ -182,19 +182,46 @@ private:
 
         auto& localProblem = *localProblems_[i];
 
-        VectorType v;
+        VectorType v, v1, v2, v3;
         localProblem.restrict(*(this->x_), v);
 
-        VectorType r;
-        localProblem.restrict(*(this->rhs_), r);
-        Dune::MatrixVector::subtractProduct(r, localProblem.getLocalMat(), v);
-        localProblem.updateLocalRhs(r);
+        v1.resize(localProblem.getLocalMat().N());
+        localProblem.getLocalMat().mv(v, v1);
+        v3.resize(matrix_.N());
+        matrix_.mv(*(this->x_), v3);
+        localProblem.restrict(v3, v2);
+
+        //print(v1, "v1");
+        //print(v2, "v2");
+
+        VectorType r1;
+        localProblem.restrict(*(this->rhs_), r1);
+        Dune::MatrixVector::subtractProduct(r1, localProblem.getLocalMat(), v);
+
+        /*localProblem.updateLocalRhs(r);
 
         VectorType x;
         localProblem.solve(x);
         v += x;
 
-        localProblem.updateVector(v, *(this->x_));
+        localProblem.updateVector(v, *(this->x_));*/
+
+        VectorType localR;
+        VectorType r = *(this->rhs_);
+        Dune::MatrixVector::subtractProduct(r, matrix_, *(this->x_));
+        localProblem.restrict(r, localR);
+
+       // print(r1, "r1");
+       // print(localR, "localR");
+
+        localProblem.updateLocalRhs(localR);
+
+        VectorType x;
+        localProblem.solve(x);
+
+        //VectorType v;
+        localProblem.prolong(x, v);
+        *(this->x_) += v;
     }
 };
 

src/cantorfaultnetworks/cantorfaultnetwork.cc

@@ -60,6 +60,8 @@
 
 #include <dune/faultnetworks/faultfactories/cantorfaultfactory.hh>
 
+#include <dune/solvers/solvers/loopsolver.hh>
+
 int main(int argc, char** argv) { try
 {
     MPIHelper::instance(argc, argv);
@@ -86,7 +88,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)
@@ -110,7 +112,6 @@ int main(int argc, char** argv) { try
     const int minCantorResolution = problemParameters.get<int>("minCantorResolution");
 
     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");
@@ -331,13 +332,23 @@ int main(int argc, char** argv) { try
     DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis);
 
     // solve
-    OscCGSolver<MatrixType, VectorType, DGBasis>
+    /*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::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN) */
+
+
+    //solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineResolution) + ".vec";
+
 
+    VectorType x = initialX;
+    preconditioner.setProblem(fineGlobalAssembler.matrix(), x, rhs);
 
-    solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineResolution) + ".vec";
+    Dune::Solvers::LoopSolver<VectorType>
+    solver (preconditioner, maxIterations, solverTolerance,
+                   fineEnergyNorm,
+                   Solver::FULL,
+                   true, &fineSol);
 
     solver.check();
     solver.preprocess();

src/cantorfaultnetworks/cantorfaultnetwork.parset

 path        = ../data/
 resultPath  = ../cantorfaultnetworks/results/
 
+[preconditioner]
+patch = SUPPORT            # CELL , SUPPORT
+mode = MULTIPLICATIVE         # ADDITIVE , MULTIPLICATIVE
+multDirection = SYMMETRIC    # SYMMETRIC , FORWARD , BACKWARD
+patchDepth = 1
 
 ###########################################
 
@@ -10,88 +15,10 @@ oscDataFile         = oscDataLaplace32.mat
 # level resolution in 2^(-...)
 coarseResolution = 0
 fineResolution = 5
-exactResolution = 10
-minCantorResolution = 0
-
-penaltyFactor = 1
-patchDepth = 1
-
-# cg solver
-nestedIteration = 0
-tol = 1e-12
-maxIterations = 8
-
-###########################################
-
-[problem1]
-oscDataFile         = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 0
-fineResolution = 6
-exactResolution = 10
+exactResolution = 8
 minCantorResolution = 0
 
 penaltyFactor = 1
-patchDepth = 1
-
-# cg solver
-nestedIteration = 0
-tol = 1e-12
-maxIterations = 8
-
-###########################################
-
-[problem2]
-oscDataFile         = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 0
-fineResolution = 7
-exactResolution = 10
-minCantorResolution = 0
-
-penaltyFactor = 1
-patchDepth = 1
-
-# cg solver
-nestedIteration = 0
-tol = 1e-12
-maxIterations = 8
-
-###########################################
-
-[problem3]
-oscDataFile         = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 0
-fineResolution = 8
-exactResolution = 10
-minCantorResolution = 0
-
-penaltyFactor = 1
-patchDepth = 1
-
-# cg solver
-nestedIteration = 0
-tol = 1e-12
-maxIterations = 8
-
-###########################################
-
-
-[problem4]
-oscDataFile         = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 0
-fineResolution = 9
-exactResolution = 10
-minCantorResolution = 0
-
-penaltyFactor = 1
-patchDepth = 1
 
 # cg solver
 nestedIteration = 0

src/cantorfaultnetworks/sparsecantorfaultnetwork.cc

@@ -59,6 +59,9 @@
 #include <dune/faultnetworks/dgmgtransfer.hh>
 
 #include <dune/faultnetworks/faultfactories/sparsecantorfaultfactory.hh>
+#include <dune/faultnetworks/faultfactories/rockfaultfactory.hh>
+
+#include <dune/solvers/solvers/loopsolver.hh>
 
 int main(int argc, char** argv) { try
 {
@@ -332,6 +335,15 @@ int main(int argc, char** argv) { try
                     maxIterations, solverTolerance, &exactEnergyNorm,
                     Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)
 
+  /*  VectorType x = initialX;
+    preconditioner.setProblem(fineGlobalAssembler.matrix(), x, rhs);
+
+    Dune::Solvers::LoopSolver<VectorType>
+    solver (preconditioner, maxIterations, solverTolerance,
+                   fineEnergyNorm,
+                   Solver::FULL,
+                   true, &fineSol);*/
+
 
     //solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec";
 

src/cantorfaultnetworks/sparsecantorfaultnetwork.parset

@@ -2,7 +2,7 @@ path        = ../data/
 resultPath  = ../cantorfaultnetworks/results/sparse/
 
 [preconditioner]
-patch = CELL            # CELL , SUPPORT
+patch = SUPPORT            # CELL , SUPPORT
 mode = MULTIPLICATIVE         # ADDITIVE , MULTIPLICATIVE
 multDirection = SYMMETRIC    # SYMMETRIC , FORWARD , BACKWARD
 patchDepth = 1
@@ -13,9 +13,9 @@ patchDepth = 1
 oscDataFile         = oscDataLaplace4.mat
 
 # level resolution in 2^(-...)
-coarseCantorLevel = 2
-fineCantorLevel = 3
-maxCantorLevel = 4
+coarseCantorLevel = 1
+fineCantorLevel = 2
+maxCantorLevel = 3
 
 penaltyFactor = 1