From 94d645d067ca81c40c02817ec36657f247183419 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 13 Nov 2014 21:38:10 +0000
Subject: [PATCH] Have an extra multigrid level when using a second- or
 third-order FE space

Previously, the number of multigrid transfer operators was always equal to
the number of grid levels minus one.  That worked out correctly when doing
first-order FE spaces.  However, for higher-order spaces, the up-most transfer
operator transfers from the higher-order space to a first-order space
_on the same grid level_.  But since the total number of transfer operators
was grid levels minus one, this implied that the coarsest grid leve
was never used, and the coarse grid solver would operate on the second-coarsest
grid.

This patch adds the extra transfer operator in this case.  That means we now
call the coarse grid solver on smaller problems, which makes a noticeable
run-time difference.

[[Imported from SVN: r9964]]
---
 dune/gfe/riemanniantrsolver.cc | 22 ++++++++++++----------
 1 file changed, 12 insertions(+), 10 deletions(-)

diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 6b3854d0..5cdb5d48 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -181,10 +181,10 @@ setup(const GridType& grid,
     for (size_t k=0; k<mmgStep->mgTransfer_.size(); k++)
         delete(mmgStep->mgTransfer_[k]);
 
-    mmgStep->mgTransfer_.resize(numLevels-1);
-
 #if defined THIRD_ORDER || defined SECOND_ORDER
-    if (numLevels>1) {
+    mmgStep->mgTransfer_.resize(numLevels);
+
+    {
         P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
 
         PKtoP1MGTransfer<CorrectionType>* topTransferOp = new PKtoP1MGTransfer<CorrectionType>;
@@ -212,19 +212,20 @@ setup(const GridType& grid,
         dynamic_cast<PKtoP1MGTransfer<CorrectionType>*>(mmgStep->mgTransfer_.back())->setMatrix(topTransferOperator);
 
         for (int i=0; i<mmgStep->mgTransfer_.size()-1; i++){
+
           // Construct the local multigrid transfer matrix
           TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
-          newTransferOp->setup(*grid_,i+1,i+2);
+          newTransferOp->setup(*grid_,i,i+1);
 
           // If we are on more than 1 processors, join all local transfer matrices on rank 0,
           // and construct a single global transfer operator there.
           typedef Dune::GlobalP1Mapper<typename GridType::LevelGridView> GlobalLevelP1Mapper;
-          GlobalLevelP1Mapper fineGUIndex(grid_->levelGridView(i+2));
-          GlobalLevelP1Mapper coarseGUIndex(grid_->levelGridView(i+1));
+          GlobalLevelP1Mapper fineGUIndex(grid_->levelGridView(i+1));
+          GlobalLevelP1Mapper coarseGUIndex(grid_->levelGridView(i));
 
           typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LevelGridView, Dune::MCMGVertexLayout> LevelLocalMapper;
-          LevelLocalMapper fineLevelLocalMapper(grid_->levelGridView(i+2));
-          LevelLocalMapper coarseLevelLocalMapper(grid_->levelGridView(i+1));
+          LevelLocalMapper fineLevelLocalMapper(grid_->levelGridView(i+1));
+          LevelLocalMapper coarseLevelLocalMapper(grid_->levelGridView(i));
 
           typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
           MatrixCommunicator<GlobalLevelP1Mapper,
@@ -232,16 +233,17 @@ setup(const GridType& grid,
                              typename GridType::LevelGridView,
                              TransferOperatorType,
                              LevelLocalMapper,
-                             LevelLocalMapper> matrixComm(fineGUIndex, coarseGUIndex, grid_->levelGridView(i+2), grid_->levelGridView(i+1), fineLevelLocalMapper, coarseLevelLocalMapper, 0);
+                             LevelLocalMapper> matrixComm(fineGUIndex, coarseGUIndex, grid_->levelGridView(i+1), grid_->levelGridView(i), fineLevelLocalMapper, coarseLevelLocalMapper, 0);
 
           mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>;
           Dune::shared_ptr<TransferOperatorType> transferOperatorMatrix = Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(newTransferOp->getMatrix()));
           dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType>*>(mmgStep->mgTransfer_[i])->setMatrix(transferOperatorMatrix);
         }
-
     }
 
 #else
+    mmgStep->mgTransfer_.resize(numLevels-1);
+
     for (size_t i=0; i<mmgStep->mgTransfer_.size(); i++){
 
         // Construct the local multigrid transfer matrix
-- 
GitLab