.

dune/faultnetworks/localproblem.hh

@@ -145,6 +145,10 @@ public:
         }
     }
 
+    void updateLocalRhs(const RangeType& rhs){
+        localRhs = rhs;
+    }
+
     void updateRhsAndBoundary(const RangeType& rhs, const RangeType& boundary) {
         for (size_t i=0; i<localRhs.size(); i++) {
             if (boundaryDofs_[i][0]) {

dune/faultnetworks/patchfactories/supportpatchfactory.hh

@@ -12,6 +12,8 @@
 
 #include <dune/faultnetworks/hierarchicleveliterator.hh>
 
+#include <dune/faultnetworks/utils/debugutils.hh>
+
 template <class BasisType>
 class SupportPatchFactory
 {
@@ -76,33 +78,68 @@ class SupportPatchFactory
 	    std::set<int> localDofs;
 	    std::set<int> localBoundaryDofs;
 
-	    addVertexSupport(centerVertexIdx_, coarseElements, visited, vertexVisited);
+        //std::set<int> centerVertices;
+        //centerVertices.insert(centerVertexIdx_);
+
+        std::set<int> coarseBoundaryElements;
+
+        addVertexSupport(centerVertexIdx_, coarseElements, visited, vertexVisited);
 	    
+        /*for (size_t i=0; i<coarseElements.size(); ++i) {
+            const auto elem = coarseElements[i];
+
+            for (const auto& is : intersections(coarseGridView, elem)) {
+                if (!is.neighbor() or is.boundary())
+                    continue;
+
+                if (visited[coarseMapper.index(is.outside())][0])
+                    continue;
+
+                const auto& localCoefficients = coarseBasis_.getLocalFiniteElement(is.inside()).localCoefficients();
+
+                for (size_t i=0; i<localCoefficients.size(); i++) {
+                    unsigned int entity = localCoefficients.localKey(i).subEntity();
+                    unsigned int codim  = localCoefficients.localKey(i).codim();
+
+                    auto idx = coarseBasis_.index(is.inside(), i);
+                    auto as = containsInsideSubentity(elem, is, entity, codim);
+                    auto asd = centerVertices.count(idx);
+
+                    if (containsInsideSubentity(elem, is, entity, codim) and centerVertices.count(coarseBasis_.index(is.inside(), i))) {
+                        coarseBoundaryElements.insert(coarseMapper.index(is.outside()));
+                        break;
+                    }
+                }
+            }
+        }*/
+
 	    // construct coarse patch
 	    for (int depth=1; depth<patchDepth_; depth++) {
-		std::vector<Element> coarseElementNeighbors;
+            std::vector<Element> coarseElementNeighbors;
 		
-		for (size_t i=0; i<coarseElements.size(); ++i) {
-		    const Element& elem = coarseElements[i];     
-		    const LFE& coarseLFE = coarseBasis_.getLocalFiniteElement(elem);
+            for (size_t i=0; i<coarseElements.size(); ++i) {
+                const Element& elem = coarseElements[i];
+                const LFE& coarseLFE = coarseBasis_.getLocalFiniteElement(elem);
 	      
-		    for (size_t j=0; j<coarseLFE.localBasis().size(); ++j) {
-            int dofIndex = coarseBasis_.indexInGridView(elem, j);
+                for (size_t j=0; j<coarseLFE.localBasis().size(); ++j) {
+                    int dofIndex = coarseBasis_.indexInGridView(elem, j);
 			
-			if (!vertexVisited[dofIndex][0]) {
-			    addVertexSupport(dofIndex, coarseElementNeighbors, visited, vertexVisited);
-			}
-		    }
-		}	        
+                    if (!vertexVisited[dofIndex][0]) {
+                        addVertexSupport(dofIndex, coarseElementNeighbors, visited, vertexVisited);
+                    }
+                }
+            }
 		    
-		coarseElements.insert(coarseElements.begin(), coarseElementNeighbors.begin(), coarseElementNeighbors.end());
+            coarseElements.insert(coarseElements.begin(), coarseElementNeighbors.begin(), coarseElementNeighbors.end());
 	    }
 	    
 
+        //print(coarseBoundaryElements, "coarseBoundaryElements");
+
 	    // construct fine patch
 	    for (size_t i=0; i<coarseElements.size(); ++i) {
 		const Element& elem = coarseElements[i];
-		addFinePatchElements(elem, visited, localDofs, localBoundaryDofs);
+        addFinePatchElements(elem, visited, localDofs, localBoundaryDofs, coarseBoundaryElements);
 	    }
 
 	 
@@ -141,13 +178,13 @@ class SupportPatchFactory
 private:
 
 
-    void print(const Dune::BitSetVector<1>& x, std::string message){
+  /*  void print(const Dune::BitSetVector<1>& x, std::string message){
        std::cout << message << std::endl;
        for (size_t i=0; i<x.size(); i++) {
            std::cout << x[i][0] << " ";
        }
        std::cout << std::endl << std::endl;
-   }
+   }*/
 
 	void addVertexSupport(const int vertexIdx, std::vector<Element>& coarseElements, Dune::BitSetVector<1>& visited, Dune::BitSetVector<1>& vertexVisited) {
 	    const std::vector<Element>& vertexInElement = vertexInElements_[vertexIdx];
@@ -169,18 +206,19 @@ private:
 	}
       
       
-	void addFinePatchElements(const Element& coarseElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
+    void addFinePatchElements(const Element& coarseElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs, const std::set<int>& coarseBoundaryElements) {
 	    HierarchicLevelIteratorType endHierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
 	    for (HierarchicLevelIteratorType hierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
 
 		const Element& fineElem = *hierIt;
 		fineRegionElements.push_back(fineElem);
 
-		addLocalDofs(coarseElem, fineElem, inCoarsePatch, localDofs, localBoundaryDofs);
+        addLocalDofs(coarseElem, fineElem, inCoarsePatch, localDofs, localBoundaryDofs, coarseBoundaryElements);
 	    }
 	}	  
 	
-	void addLocalDofs(const Element& coarseElem, const Element& fineElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
+    void addLocalDofs(const Element& coarseElem, const Element& fineElem, const Dune::BitSetVector<1>& inCoarsePatch,
+                      std::set<int>& localDofs, std::set<int>& localBoundaryDofs, const std::set<int>& coarseBoundaryElements) {
 	    const LFE& fineLFE = fineBasis_.getLocalFiniteElement(fineElem);
 	    
         /*
@@ -210,8 +248,9 @@ private:
 			outsideFather = outsideFather.father();
 		    }
 		    
-		    if (!inCoarsePatch[coarseMapper.index(outsideFather)][0]) {
-			isInnerBoundary = true;
+            auto outsideFatherIdx = coarseMapper.index(outsideFather);
+            if (!inCoarsePatch[outsideFatherIdx][0] and !coarseBoundaryElements.count(outsideFatherIdx)) {
+                isInnerBoundary = true;
 		    }
         }
 		
@@ -219,9 +258,9 @@ private:
 		    typedef typename LFE::Traits::LocalCoefficientsType LocalCoefficients;
 		    const LocalCoefficients* localCoefficients = &fineBasis_.getLocalFiniteElement(intersection.inside()).localCoefficients();
 
-		    for (size_t i=0; i<localCoefficients->size(); i++) {
-			unsigned int entity = localCoefficients->localKey(i).subEntity();
-			unsigned int codim  = localCoefficients->localKey(i).codim();
+            for (size_t i=0; i<localCoefficients->size(); i++) {
+            unsigned int entity = localCoefficients->localKey(i).subEntity();
+            unsigned int codim  = localCoefficients->localKey(i).codim();
 
 			if (containsInsideSubentity(fineElem, intersection, entity, codim))
 			    localBoundaryDofs.insert(fineBasis_.index(intersection.inside(), i));

dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh

@@ -186,7 +186,9 @@ private:
         localProblem.restrict(*(this->x_), v);
 
         VectorType r;
-        localProblem.localResidual(v, r);
+        localProblem.restrict(*(this->rhs_), r);
+        Dune::MatrixVector::subtractProduct(r, localProblem.getLocalMat(), v);
+        localProblem.updateLocalRhs(r);
 
         VectorType x;
         localProblem.solve(x);

dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh

@@ -6,6 +6,8 @@
 #include <dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh>
 #include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
 
+#include <dune/faultnetworks/utils/debugutils.hh>
+
 template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
 class SupportPatchPreconditioner : public LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType> {
 private:
@@ -31,7 +33,7 @@ public:
         std::vector<std::vector<Element>> vertexInElements;
 
         const GridView& patchLevelGridView = this->patchLevelBasis_.getGridView();
-        vertexInElements.resize(patchLevelGridView.size(dim));
+        vertexInElements.resize(patchLevelGridView.size(dim)); //(this->patchLevelBasis_.size());
         for (size_t i=0; i<vertexInElements.size(); i++) {
             vertexInElements[i].resize(0);
         }
@@ -41,13 +43,15 @@ public:
             const auto& coarseFE = this->patchLevelBasis_.getLocalFiniteElement(elem);
 
             for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
-                int dofIndex = this->patchLevelBasis_.indexInGridView(elem, i);
+                int dofIndex = this->patchLevelBasis_.indexInGridView(elem, i); //index
                 vertexInElements[dofIndex].push_back(elem);
             }
         }
 
         std::cout << "SupportPatchPreconditioner::build() level: " << this->level_ << std::endl;
 
+        //printBasisDofLocation(this->patchLevelBasis_);
+
         this->localProblems_.resize(vertexInElements.size());
         const int patchLevel = this->patchLevelBasis_.faultNetwork().level();
         for (size_t i=0; i<vertexInElements.size(); i++) {

dune/faultnetworks/solvers/osccgsolver.cc

@@ -232,10 +232,6 @@ void OscCGSolver<MatrixType, VectorType, DGBasisType>::solve()
         std::cout << "--------------------\n";
 
 
-
-	if (this->useRelativeError_)
-            std::cout << "Norm of final discretization error: " << (discretizationError*normOfInitialDiscretizationError) << std::endl;  
-        else
-            std::cout << "Norm of final discretization error: " << discretizationError << std::endl;  
+       std::cout << "Norm of final discretization error: " << discretizationError << std::endl;
     }
 }

src/cantorfaultnetworks/sparsecantorfaultnetwork.cc

@@ -333,7 +333,7 @@ int main(int argc, char** argv) { try
                     Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)
 
 
-    solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec";
+    //solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec";
 
     solver.check();
     solver.preprocess();

src/cantorfaultnetworks/sparsecantorfaultnetwork.parset

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