.

dune/faultnetworks/patchfactories/cellpatchfactory.hh

@@ -23,9 +23,7 @@ class CellPatchFactory
 	typedef typename GridType::ctype ctype;
 	static const int dim = GridType::dimension;
 	typedef typename GridType::template Codim<0>::Entity Element;
-	typedef typename GridType::template Codim<dim>::Entity Vertex;
 
-	typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
 	typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,  Dune::MCMGElementLayout > ElementMapper;
 	
     private:
@@ -42,13 +40,13 @@ class CellPatchFactory
     ElementMapper mapper;
 
 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;
-    }
+    }*/
 
     bool containsInsideSubentity(const Element& elem, const typename GridType::LevelIntersection& intersection, int subEntity, int codim) {
         return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);

dune/faultnetworks/patchfactories/supportpatchfactory.hh

 #ifndef SUPPORT_PATCH_FACTORY_HH
 #define SUPPORT_PATCH_FACTORY_HH
 
-#include<queue>
+#include <queue>
+#include <set>
 
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/fvector.hh>
@@ -227,27 +228,6 @@ private:
 		    } 
 		}
 	    }	
-	    
-        /* add interior coarse level vertices to boundary dofs
-	    for(int i=0; i<coarseGeometry.corners(); ++i) {
-		const auto& local = fineGeometry.local(coarseGeometry.corner(i));
-		    
-		if (!ref.checkInside(local))
-		    continue;
-		
-		const int localBasisSize = fineLFE.localBasis().size(); 
-		std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
-		fineLFE.localBasis().evaluateFunction(local, localBasisRep);
-		
-		for(int k=0; k<localBasisSize; ++k) {
-		    if (localBasisRep[k]==1) {
-			int dofIndex = fineBasis_.index(fineElem, k);
-			localBoundaryDofs.insert(dofIndex);
-			break;
-		    }	
-		}
-        }   */
-
 	}
 };
 #endif

dune/faultnetworks/preconditioners/CMakeLists.txt

@@ -2,12 +2,16 @@
 
 add_custom_target(faultnetworks_preconditioners_src SOURCES
   levelglobalpreconditioner.hh
+  cellpatchpreconditioner.hh
+  supportpatchpreconditioner.hh
   levelpatchpreconditioner.hh
   multilevelpatchpreconditioner.hh
 )
 
 install(FILES
   levelglobalpreconditioner.hh
+  cellpatchpreconditioner.hh
+  supportpatchpreconditioner.hh
   levelpatchpreconditioner.hh
   multilevelpatchpreconditioner.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/preconditioners)

dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh

+#ifndef CELL_PATCH_PRECONDITIONER_HH
+#define CELL_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh>
+#include <dune/faultnetworks/patchfactories/cellpatchfactory.hh>
+
+template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+class CellPatchPreconditioner : public LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType> {
+private:
+    using Base = LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+    using GridView = typename Base::GridView;
+    using Mode = typename Base::Mode;
+
+public:
+    CellPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+                             const LocalAssembler& localAssembler,
+                             const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+                             const Mode mode = Mode::additive) :
+          Base(levelInterfaceNetwork, localAssembler, localInterfaceAssemblers, mode) {}
+
+    void build() {
+        CellPatchFactory<BasisType> patchFactory(this->basis_, this->level_);
+        const auto& localToGlobal = patchFactory.getLocalToGlobal();
+        const auto& boundaryDofs = patchFactory.getBoundaryDofs();
+
+        VectorType rhs;
+        rhs.resize(this->matrix_.N());
+        rhs = 0;
+
+        std::cout << "CellPatchPreconditioner::build() level: " << this->level_ << std::endl;
+
+        size_t cellCount = localToGlobal.size();
+        this->localProblems_.resize(cellCount);
+        for (size_t i=0; i<cellCount; i++) {
+            this->localProblems_[i] = std::make_shared<decltype(*(this->localProblems_[i]))>(this->matrix_, rhs, localToGlobal[i], boundaryDofs[i]);
+        }
+    }
+};
+
+#endif
+

dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh

@@ -10,7 +10,6 @@
 #include <dune/solvers/iterationsteps/lineariterationstep.hh>
 
 #include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
-#include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
 #include <dune/faultnetworks/localproblem.hh>
 #include <dune/faultnetworks/levelinterfacenetwork.hh>
 #include <dune/faultnetworks/utils/debugutils.hh>
@@ -18,20 +17,19 @@
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functiontools/boundarydofs.hh>
 
-
 template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
 class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
 
 public:
-    enum Mode {additive, multiplicative};
+    enum Mode {ADDITIVE, MULTIPLICATIVE};
     enum BoundaryMode {homogeneous, fromIterate};
+    enum Direction {FORWARD, BACKWARD, SYMMETRIC};
 
 private:
     typedef typename BasisType::GridView GridView;
     typedef typename GridView::Grid GridType;
 
     const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
-    const BasisType& patchLevelBasis_;
     const LocalAssembler& localAssembler_;
     const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
     const Mode mode_;
@@ -46,30 +44,8 @@ private:
     MatrixType matrix_;
     std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
 
-public:
-
-    // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
-    LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
-                             const BasisType& patchLevelBasis,
-                             const LocalAssembler& localAssembler,
-                             const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
-                             const Mode mode = LevelPatchPreconditioner::Mode::additive) :
-          levelInterfaceNetwork_(levelInterfaceNetwork),
-          patchLevelBasis_(patchLevelBasis),
-          localAssembler_(localAssembler),
-          localInterfaceAssemblers_(localInterfaceAssemblers),
-          mode_(mode),
-          grid_(levelInterfaceNetwork_.grid()),
-          level_(levelInterfaceNetwork_.level()),
-          basis_(levelInterfaceNetwork_)
-    {
-        setPatchDepth();
-        setBoundaryMode();
-
-        assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
-    }
-
-    void build() {
+private:
+    void setup() {
         // assemble stiffness matrix for entire level including all faults
         GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
         globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
@@ -96,70 +72,31 @@ public:
             }
             row[i] = 1;
         }
+    }
 
-        // init vertexInElements
-        const int dim = GridType::dimension;
-        typedef typename GridType::template Codim<0>::Entity EntityType;
-        std::vector<std::vector<EntityType>> vertexInElements;
-
-        const GridView& patchLevelGridView = patchLevelBasis_.getGridView();
-        vertexInElements.resize(patchLevelGridView.size(dim));
-        for (size_t i=0; i<vertexInElements.size(); i++) {
-            vertexInElements[i].resize(0);
-        }
-
-        typedef typename BasisType::LocalFiniteElement FE;
-        typedef  typename GridView::template Codim <0>::Iterator  ElementLevelIterator;
-        ElementLevelIterator endElemIt = patchLevelGridView.template end <0>();
-        for (ElementLevelIterator  it = patchLevelGridView. template begin <0>(); it!=endElemIt; ++it) {
-
-            // compute coarseGrid vertexInElements
-            const FE& coarseFE = patchLevelBasis_.getLocalFiniteElement(*it);
-
-            for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
-                int dofIndex = patchLevelBasis_.indexInGridView(*it, i);
-                vertexInElements[dofIndex].push_back(*it);
-            }
-        }
-
-        localProblems_.resize(vertexInElements.size());
-
-        std::cout << std::endl;
-        std::cout << "---------------------------------------------" << std::endl;
-        std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;
-
-        // init local fine level corrections
-        Dune::Timer timer;
-        timer.reset();
-        timer.start();
-
-        const int patchLevel = patchLevelBasis_.faultNetwork().level();
-        for (size_t i=0; i<vertexInElements.size(); i++) {
-            //std::cout << i << std::endl;
-            //std::cout << "---------------" << std::endl;
-
-            SupportPatchFactory<BasisType> patchFactory(patchLevelBasis_, basis_, patchLevel, level_, vertexInElements, i, patchDepth_);
-            std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
-            Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
-            VectorType rhs;
-            rhs.resize(matrix_.N());
-            rhs = 0;
-
-            //print(localToGlobal, "localToGlobal: ");
-            //print(boundaryDofs, "boundaryDofs: ");
+public:
+    // has to setup localProblems_
+    virtual void build();
 
-            localProblems_[i] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs, localToGlobal, boundaryDofs);
+public:
+    LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+                             const LocalAssembler& localAssembler,
+                             const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+                             const Mode mode = LevelPatchPreconditioner::Mode::additive) :
+          levelInterfaceNetwork_(levelInterfaceNetwork),
+          localAssembler_(localAssembler),
+          localInterfaceAssemblers_(localInterfaceAssemblers),
+          mode_(mode),
+          grid_(levelInterfaceNetwork_.grid()),
+          level_(levelInterfaceNetwork_.level()),
+          basis_(levelInterfaceNetwork_)
+    {
+        setPatchDepth();
+        setBoundaryMode();
 
-            /*
-            if ((i+1) % 10 == 0) {
-                std::cout << '\r' << std::floor((i+1.0)/patchLevelBasis_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*patchLevelBasis_.size() << " seconds." << std::flush;
-            }*/
-        }
+        assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
 
-        std::cout << std::endl;
-        std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
-        std::cout << "---------------------------------------------" << std::endl;
-        timer.stop();
+        setup();
     }
 
     void setPatchDepth(const size_t patchDepth = 0) {
@@ -183,13 +120,30 @@ public:
         }
     }
 
-    virtual void iterate() {
-        if (mode_ == additive)
+    virtual void iterate(Direction dir = SYMMETRIC) {
+        if (mode_ == ADDITIVE)
             iterateAdd();
         else
-            iterateMult();
+            iterateMult(dir);
+    }
+
+    const BasisType& basis() const {
+        return basis_;
+    }
+
+    const GridView& gridView() const {
+        return basis_.getGridView();
+    }
+
+    const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
+        return levelInterfaceNetwork_;
+    }
+
+    size_t size() const {
+        return localProblems_.size();
     }
 
+private:
     void iterateAdd() {
         *(this->x_) = 0;	
 
@@ -206,54 +160,38 @@ public:
         }
     }
 
-    void iterateMult() {
+    void iterateMult(Direction dir) {
+        if (dir != BACKWARD) {
+            for (size_t i=0; i<localProblems_.size(); i++)
+                iterateStep(i);
+        }
+
+        if (dir != Direction::FORWARD)
+            for (size_t i=localProblems_.size()-1; i>=0 && i<localProblems_.size(); i--)
+                iterateStep(i);
+    }
+
+    void iterateStep(size_t i) {
         *(this->x_) = 0;
 
-        VectorType it, x;
-        for (size_t i=0; i<localProblems_.size(); i++) {
-            VectorType updatedRhs(*(this->rhs_));
-            matrix_.mmv(*(this->x_), updatedRhs);
+        auto rhs = *(this->rhs_);
 
-            //step(i);
-            //print(updatedRhs, "localRhs: ");
-            //writeToVTK(basis_, updatedRhs, "/storage/mi/podlesny/data/faultnetworks/rhs/local", "exactvertexdata_step_"+std::to_string(i));
+        VectorType it, x;
 
             if (boundaryMode_ == BoundaryMode::homogeneous)
-                localProblems_[i]->updateRhs(updatedRhs);
+                localProblems_[i]->updateRhs(rhs);
             else
-                localProblems_[i]->updateRhsAndBoundary(updatedRhs, *(this->x_));
+                localProblems_[i]->updateRhsAndBoundary(rhs, *(this->x_));
 
             localProblems_[i]->solve(it);
             localProblems_[i]->prolong(it, x);
 
+            VectorType residual;
+            localProblems_[i]->localResidual(it, residual);
+            localProblems_[i]->updateVector(residual, rhs);
 
-            //print(it, "it: ");
-            //print(x, "x: ");
-
-            //writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/sol/local", "exactvertexdata_step_"+std::to_string(i));
-
-            /*if (i==5) {
-                writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
-            }*/
 
             *(this->x_) += x;
-        }
-    }
-
-    const BasisType& basis() const {
-        return basis_;
-    }
-
-    const GridView& gridView() const {
-        return basis_.getGridView();
-    }
-
-    const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
-        return levelInterfaceNetwork_;
-    }
-
-    size_t size() const {
-        return localProblems_.size();
     }
 };
 

dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh

+#ifndef SUPPORT_PATCH_PRECONDITIONER_HH
+#define SUPPORT_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh>
+#include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
+
+template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+class SupportPatchPreconditioner : public LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType> {
+private:
+    using Base = LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+    using GridView = typename Base::GridView;
+    using GridType = typename Base::GridType;
+    using Mode = typename Base::Mode;
+
+    const BasisType& patchLevelBasis_;
+
+public:
+    // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
+    SupportPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+                             const BasisType& patchLevelBasis,
+                             const LocalAssembler& localAssembler,
+                             const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+                             const Mode mode = Mode::additive) :
+          Base(levelInterfaceNetwork, localAssembler, localInterfaceAssemblers, mode),
+          patchLevelBasis_(patchLevelBasis) {}
+
+    void build() {
+        // init vertexInElements
+        const int dim = GridType::dimension;
+        using Element = typename GridType::template Codim<0>::Entity;
+        std::vector<std::vector<Element>> vertexInElements;
+
+        const GridView& patchLevelGridView = patchLevelBasis_.getGridView();
+        vertexInElements.resize(patchLevelGridView.size(dim));
+        for (size_t i=0; i<vertexInElements.size(); i++) {
+            vertexInElements[i].resize(0);
+        }
+
+        for (const auto& elem : elements(patchLevelGridView)) {
+            // compute coarseGrid vertexInElements
+            const auto& coarseFE = patchLevelBasis_.getLocalFiniteElement(elem);
+
+            for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
+                int dofIndex = patchLevelBasis_.indexInGridView(elem, i);
+                vertexInElements[dofIndex].push_back(elem);
+            }
+        }
+
+        std::cout << "SupportPatchPreconditioner::build() level: " << this->level_ << std::endl;
+
+        this->localProblems_.resize(vertexInElements.size());
+        const int patchLevel = patchLevelBasis_.faultNetwork().level();
+        for (size_t i=0; i<vertexInElements.size(); i++) {
+            SupportPatchFactory<BasisType> patchFactory(patchLevelBasis_, this->basis_, patchLevel, this->level_, vertexInElements, i, this->patchDepth_);
+            const auto& localToGlobal = patchFactory.getLocalToGlobal();
+            const auto& boundaryDofs = patchFactory.getBoundaryDofs();
+
+            VectorType rhs;
+            rhs.resize(this->matrix_.N());
+            rhs = 0;
+
+            this->localProblems_[i] = std::make_shared<decltype(*(this->localProblems_[i]))>(this->matrix_, rhs, localToGlobal, boundaryDofs);
+        }
+    }
+};
+
+#endif
+

src/cantorfaultnetworks/sparsecantorfaultnetwork.cc

@@ -302,7 +302,7 @@ int main(int argc, char** argv) { try
 
     typedef MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelPatchPreconditioner;
     typedef typename MultilevelPatchPreconditioner::LevelPatchPreconditionerType LevelPatchPreconditioner;
-    MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::multiplicative;
+    MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::additive;
 
     MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
     for (size_t i=0; i<preconditioner.size()-1; i++) {