Allow to set matrix hierarchy manually in MultiGridStep

This is helpful when using AMG, where the matrix hierarchy
is computed while coarsening. Another application are
problems, where the coarse matrix is a preconditioner
that does not coincide with the Galerkin-restriction
of the fine matrix, e.g. for DG methods.

dune/solvers/iterationsteps/multigridstep.cc

@@ -73,17 +73,14 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
 
     xHierarchy_.resize(numLevels());
     rhsHierarchy_.resize(numLevels());
-    matrixHierarchy_.resize(numLevels());
     ignoreNodesHierarchy_.resize(numLevels());
 
-    for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
+    for (int i=0; i<int(xHierarchy_.size())-1; i++)
     {
-        matrixHierarchy_[i].reset();
         xHierarchy_[i].reset();
         ignoreNodesHierarchy_[i] = NULL;
     }
 
-    matrixHierarchy_.back() = this->mat_;
     xHierarchy_.back()      = Dune::stackobject_to_shared_ptr(*(this->x_));
     rhsHierarchy_.back()    = *(this->rhs_);
 
@@ -117,22 +114,35 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
     // /////////////////////////////////////////////////////
     //   Assemble the complete hierarchy of matrices
     // /////////////////////////////////////////////////////
-    for (int i=this->numLevels()-2; i>=0; i--)
+    if (not matrixHierarchyManuallySet_)
     {
-        this->matrixHierarchy_[i] = std::make_shared<MatrixType>();
-        this->xHierarchy_[i] = std::make_shared<VectorType>();
+        matrixHierarchy_.resize(numLevels());
+        for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
+            matrixHierarchy_[i].reset();
+        matrixHierarchy_.back() = this->mat_;
+        for (int i=this->numLevels()-2; i>=0; i--)
+        {
+            this->matrixHierarchy_[i] = std::make_shared<MatrixType>();
 
-        // Compute which entries are present in the (sparse) coarse grid stiffness
-        // matrices.
-        this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+            // Compute which entries are present in the (sparse) coarse grid stiffness
+            // matrices.
+            this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+
+            // setup matrix
+            this->mgTransfer_[i]->galerkinRestrict(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+        }
+        matrixHierarchyManuallySet_ = false;
+    }
 
-        // setup matrix
-        this->mgTransfer_[i]->galerkinRestrict(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+    for (int i=this->numLevels()-2; i>=0; i--)
+    {
+        this->xHierarchy_[i] = std::make_shared<VectorType>();
 
         // Set solution vector sizes for the lower levels
         MatrixVector::resize(*(this->xHierarchy_[i]),*this->matrixHierarchy_[i]);
     }
 
+
     // /////////////////////////////////////////////////////
     //   Setup dirichlet bitfields
     // /////////////////////////////////////////////////////

dune/solvers/iterationsteps/multigridstep.hh

@@ -208,6 +208,24 @@ namespace Dune {
             basesolver_ = wrap_own_share<Solver>(std::forward<BaseSolver>(baseSolver));
         }
 
+        /**
+         * \brief Set pre-coarsened matrix hierarchy
+         *
+         * This allows to set the coarsened matrix hierarchy manually.
+         * It is assumed, that the finest matrix coincides with
+         * the one passed to setProblem(). The latter is ignored.
+         * This has to be called before preprocess().
+         */
+        template<class M>
+        void setMatrixHirarchy(const std::vector<std::shared_ptr<M>>& matrixHierarchy)
+        {
+            matrixHierarchy_.resize(matrixHierarchy.size());
+            for(auto i : Dune::range(matrixHierarchy.size()))
+                matrixHierarchy_[i] = matrixHierarchy[i];
+            matrixHierarchyManuallySet_ = true;
+        }
+
+
     protected:
         /** \brief The presmoothers, one for each level */
         std::vector<std::shared_ptr<LinearStepType> > presmoother_;
@@ -226,6 +244,8 @@ namespace Dune {
 
         //! Number of coarse corrections
         int mu_;
+
+        bool matrixHierarchyManuallySet_= false;
     public:
         //! Variable used to store the current level
         int level_;