Implement (algebraic) nested iteration

This purely algebraic variant of nested iteration
can be used to compute initial iterates.

[[Imported from SVN: r5920]]

dune/solvers/iterationsteps/obstacletnnmgstep.hh

@@ -17,6 +17,7 @@
 #include <dune/solvers/solvers/loopsolver.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
 
 
 
@@ -335,6 +336,83 @@ class ObstacleTNNMGStep
             return s;
         }
 
+
+        /**
+         * \brief Compute initial iterate by nested iteration
+         *
+         * This method performs a purely algebraic variant of nested iteration.
+         * The coarse matrix and rhs are obtained using the transfer operators.
+         * Coarse obstacles are constructed using a local max/min over the
+         * support of coarse basis functions as suggested by J. Mandel.
+         * The ignore-marks are restricted such that a coarse dof
+         * is ignored if it is associated to an ignored fine dof
+         * in the sense that the corresponding transfer operators entry is 1.
+         *
+         * On each level a fixed number of v-cycles is performed.
+         *
+         * This method can be called before the preprocess() method.
+         * It does not change the ObstacleTNNMGStep state despite
+         * updating the solution vector.
+         *
+         * \param coarseIterationSteps Number of v-cycle performed on the corse levels.
+         */
+        void nestedIteration(int coarseIterationSteps=2)
+        {
+            int maxLevel = transfer_.size();
+
+            std::vector<Matrix> coarseMatrix(maxLevel);
+            std::vector<Vector> coarseSolution(maxLevel);
+            std::vector<Vector> coarseRhs(maxLevel);
+            std::vector<BitVector> coarseIgnore(maxLevel);
+            std::vector<BoxConstraintVector> coarseObstacle(maxLevel);
+            MandelObstacleRestrictor<Vector> obstacleRestrictor;
+            BitVector critical;
+
+            critical.resize(x_->size());
+            critical.unsetAll();
+
+            transfer_[maxLevel-1]->galerkinRestrictSetOccupation(mat_, coarseMatrix[maxLevel-1]);
+            transfer_[maxLevel-1]->galerkinRestrict(mat_, coarseMatrix[maxLevel-1]);
+            transfer_[maxLevel-1]->restrict(rhs_, coarseRhs[maxLevel-1]);
+            transfer_[maxLevel-1]->restrictToFathers(*ignoreNodes_, coarseIgnore[maxLevel-1]);
+            obstacleRestrictor.restrict(obstacles_, coarseObstacle[maxLevel-1], hasObstacle_, hasObstacle_, *(transfer_[maxLevel-1]), critical);
+            for (int i = maxLevel-2; i>=0; --i)
+            {
+                transfer_[i]->galerkinRestrictSetOccupation(coarseMatrix[i+1], coarseMatrix[i]);
+                transfer_[i]->galerkinRestrict(coarseMatrix[i+1], coarseMatrix[i]);
+                transfer_[i]->restrict(coarseRhs[i+1], coarseRhs[i]);
+                transfer_[i]->restrictToFathers(coarseIgnore[i+1], coarseIgnore[i]);
+                obstacleRestrictor.restrict(coarseObstacle[i+1], coarseObstacle[i], hasObstacle_, hasObstacle_, *(transfer_[i]), critical);
+                coarseSolution[i].resize(coarseMatrix[i].N());
+            }
+
+            coarseSolution[0] = 0;
+            for (int i = 0; i<maxLevel; ++i)
+            {
+
+                TransferPointerVector coarseTransfer(i);
+                for (int j = 0; j < i; ++j)
+                    coarseTransfer[j] = transfer_[j];
+
+                ObstacleTNNMGStep coarseTNNMGStep(
+                        coarseMatrix[i],
+                        coarseSolution[i],
+                        coarseRhs[i],
+                        coarseIgnore[i],
+                        coarseObstacle[i],
+                        coarseTransfer,
+                        truncationTol_);
+                coarseTNNMGStep.preprocess();
+                for (int j = 0; j < coarseIterationSteps; ++j)
+                    coarseTNNMGStep.iterate();
+                coarseSolution[i] = coarseTNNMGStep.getSol();
+
+                if (i<maxLevel-1)
+                    transfer_[i]->prolong(coarseSolution[i], coarseSolution[i+1]);
+            }
+            transfer_[maxLevel-1]->prolong(coarseSolution[maxLevel-1], *x_);
+        }
+
     protected:
 
         // problem description given by the user