I compute the local-to-global index mapping within the MatrixCommunicator

In particular this means having to carry mappers for the local indices around.
This makes the code more verbose, but it also makes it easier to see what is
happening.  While coding this stuff I keep having doubts about whether all
the renumbering works correctly...

[[Imported from SVN: r9925]]

dune/gfe/parallel/globalp2mapper.hh

@@ -32,17 +32,16 @@ namespace Dune {
     typedef std::map<Index,Index>    IndexMap;
 
     GlobalP2Mapper(const GridView& gridView)
+    : p2Mapper_(gridView)
     {
       static_assert(GridView::dimension==2, "Only implemented for two-dimensional grids");
 
-      P2BasisMapper<GridView> p2Mapper(gridView);
-
       GlobalIndexSet<GridView> globalVertexIndex(gridView,2);
       GlobalIndexSet<GridView> globalEdgeIndex(gridView,1);
       GlobalIndexSet<GridView> globalElementIndex(gridView,0);
 
       // total number of degrees of freedom
-      nGlobalEntity_ = globalVertexIndex.nGlobalEntity() + globalEdgeIndex.nGlobalEntity() + globalElementIndex.nGlobalEntity();
+      nGlobalEntity_ = globalVertexIndex.size(2) + globalEdgeIndex.size(1) + globalElementIndex.size(0);
       nOwnedLocalEntity_ = globalVertexIndex.nOwnedLocalEntity() + globalEdgeIndex.nOwnedLocalEntity() + globalElementIndex.nOwnedLocalEntity();
 
       // Determine
@@ -56,7 +55,7 @@ namespace Dune {
 #endif
         {
           //int localIndex  = globalVertexIndex.localIndex (*it->template subEntity<2>(i));
-          int localIndex  = p2Mapper.map(*it, i, 2);
+          int localIndex  = p2Mapper_.map(*it, i, 2);
           int globalIndex = globalVertexIndex.index(*it->template subEntity<2>(i));
 
           localGlobalMap_[localIndex]  = globalIndex;
@@ -71,8 +70,8 @@ namespace Dune {
 #endif
         {
           //int localIndex  = globalEdgeIndex.localIndex (*it->template subEntity<1>(i)) + gridView.size(2);
-          int localIndex  = p2Mapper.map(*it, i, 1);
-          int globalIndex = globalEdgeIndex.index(*it->template subEntity<1>(i)) + globalVertexIndex.nGlobalEntity();
+          int localIndex  = p2Mapper_.map(*it, i, 1);
+          int globalIndex = globalEdgeIndex.index(*it->template subEntity<1>(i)) + globalVertexIndex.size(2);
 
           localGlobalMap_[localIndex]  = globalIndex;
           globalLocalMap_[globalIndex] = localIndex;
@@ -85,7 +84,7 @@ namespace Dune {
         if (it->type().isQuadrilateral())
         {
           //int localIndex  = globalEdgeIndex.localIndex (*it->template subEntity<1>(i)) + gridView.size(2);
-          int localIndex  = p2Mapper.map(*it, 0, 0);
+          int localIndex  = p2Mapper_.map(*it, 0, 0);
           int globalIndex = globalElementIndex.index(*it->template subEntity<0>(0))
                             + globalEdgeIndex.nGlobalEntity()
                             + globalVertexIndex.nGlobalEntity();
@@ -103,6 +102,13 @@ namespace Dune {
       return localGlobalMap_.find(localIndex)->second;
     }
 
+    template <class Entity>
+    Index subIndex(const Entity& entity, uint i, uint codim) const
+    {
+      int localIndex = p2Mapper_.map(entity, i, codim);
+      return localGlobalMap_.find(localIndex)->second;
+    }
+
     Index localIndex(const int& globalIndex) const {
       return globalLocalMap_.find(globalIndex)->second;
     }
@@ -117,6 +123,8 @@ namespace Dune {
       return nOwnedLocalEntity_;
     }
 
+    P2BasisMapper<GridView> p2Mapper_;
+
     IndexMap localGlobalMap_;
     IndexMap globalLocalMap_;
 

dune/gfe/parallel/matrixcommunicator.hh

@@ -9,7 +9,7 @@
 #include <dune/gfe/parallel/mpifunctions.hh>
 
 
-template<typename GUIndex, typename Communicator, typename MatrixType, typename ColGUIndex=GUIndex>
+template<typename GUIndex, typename GridView, typename MatrixType, typename LocalMapper1, typename LocalMapper2, typename ColGUIndex=GUIndex>
 class MatrixCommunicator {
 
   struct TransferMatrixTuple {
@@ -34,7 +34,7 @@ class MatrixCommunicator {
         const int i = rIt.index();
         const int j = cIt.index();
 
-        localMatrixEntries.push_back(TransferMatrixTuple(guIndex1_.index(i), guIndex2_.index(j), *cIt));
+        localMatrixEntries.push_back(TransferMatrixTuple(localToGlobal1_[i], localToGlobal2_[j], *cIt));
       }
 
     // Get number of matrix entries on each process
@@ -45,19 +45,27 @@ class MatrixCommunicator {
   }
 
 public:
-  MatrixCommunicator(const GUIndex& rowIndex, const Communicator& communicator, const int& root)
+  MatrixCommunicator(const GUIndex& rowIndex, const GridView& gridView, const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root)
   : guIndex1_(rowIndex),
     guIndex2_(rowIndex),
-    communicator_(communicator),
+    localMapper1_(localMapper1),
+    localMapper2_(localMapper2),
+    communicator_(gridView.comm()),
     root_rank(root)
-  {}
+  {
+    setLocalToGlobal(gridView);
+  }
 
-  MatrixCommunicator(const GUIndex& rowIndex, const ColGUIndex& colIndex, const Communicator& communicator, const int& root)
+  MatrixCommunicator(const GUIndex& rowIndex, const ColGUIndex& colIndex, const GridView& gridView, const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root)
   : guIndex1_(rowIndex),
     guIndex2_(colIndex),
-    communicator_(communicator),
+    localMapper1_(localMapper1),
+    localMapper2_(localMapper2),
+    communicator_(gridView.comm()),
     root_rank(root)
-  {}
+  {
+    setLocalToGlobal(gridView);
+  }
 
   MatrixType reduceAdd(const MatrixType& local)
   {
@@ -110,11 +118,42 @@ public:
   }
 
 private:
+
+  void setLocalToGlobal(const GridView& gridView)
+  {
+    localToGlobal1_.resize(localMapper1_.size());
+    localToGlobal2_.resize(localMapper2_.size());
+
+    for (auto it = gridView.template begin<0>(); it != gridView.template end<0>(); ++it)
+      for (int codim = 0; codim <= GridView::dimension; codim++)
+        for (size_t i=0; i<it->subEntities(codim); i++)
+        {
+          typename GUIndex::Index localIdx = localMapper1_.map(*it,i,codim);
+          typename GUIndex::Index globalIdx = guIndex1_.subIndex(*it,i,codim);
+          localToGlobal1_[localIdx] = globalIdx;
+
+          localIdx = localMapper2_.map(*it,i,codim);
+          globalIdx = guIndex2_.subIndex(*it,i,codim);
+          localToGlobal2_[localIdx] = globalIdx;
+        }
+
+
+  }
+
+  // Mappers for the global numbering
   const GUIndex& guIndex1_;
   const ColGUIndex& guIndex2_;
-  const Communicator& communicator_;
+
+  // Mappers for the local numbering
+  const LocalMapper1& localMapper1_;
+  const LocalMapper2& localMapper2_;
+
+  const typename GridView::CollectiveCommunication& communicator_;
   int root_rank;
 
+  std::vector<typename GUIndex::Index> localToGlobal1_;
+  std::vector<typename ColGUIndex::Index> localToGlobal2_;
+
   std::vector<TransferMatrixTuple> globalMatrixEntries;
 };
 

dune/gfe/riemanniantrsolver.cc

@@ -6,6 +6,8 @@
 
 #include <dune/istl/io.hh>
 
+#include <dune/grid/common/mcmgmapper.hh>
+
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 #include <dune/fufem/assemblers/operatorassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
@@ -126,8 +128,14 @@ setup(const GridType& grid,
     if (h1SemiNorm_)
         delete h1SemiNorm_;
 
+    LocalMapper localMapper(grid_->leafGridView());
+
+    MatrixCommunicator<GUIndex,
+                       typename GridType::LeafGridView,
+                       ScalarMatrixType,
+                       LocalMapper,
+                       LocalMapper> matrixComm(*guIndex_, grid_->leafGridView(), localMapper, localMapper, 0);
 
-    MatrixCommunicator<GUIndex, typename GridType::LeafGridView::CollectiveCommunication, ScalarMatrixType> matrixComm(*guIndex_, grid_->leafGridView().comm(), 0);
     ScalarMatrixType* A = new ScalarMatrixType(matrixComm.reduceAdd(localA));
 
     h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
@@ -186,11 +194,16 @@ setup(const GridType& grid,
         typedef Dune::GlobalIndexSet<typename GridType::LeafGridView> LeafP1GUIndex;
         LeafP1GUIndex p1Index(grid_->leafGridView(), gridDim);
 
+        typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView, Dune::MCMGVertexLayout> LeafP1LocalMapper;
+        LeafP1LocalMapper leafP1LocalMapper(grid_->leafGridView());
+
         typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
         MatrixCommunicator<GUIndex,
-                           typename GridType::LeafGridView::CollectiveCommunication,
+                           typename GridType::LeafGridView,
                            TransferOperatorType,
-                           LeafP1GUIndex> matrixComm(*guIndex_, p1Index, grid_->leafGridView().comm(), 0);
+                           LocalMapper,
+                           LeafP1LocalMapper,
+                           LeafP1GUIndex> matrixComm(*guIndex_, p1Index, grid_->leafGridView(), localMapper, leafP1LocalMapper, 0);
 
         mmgStep->mgTransfer_.back() = new PKtoP1MGTransfer<CorrectionType>;
         Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(topTransferOp->getMatrix()));
@@ -207,10 +220,16 @@ setup(const GridType& grid,
           LevelGUIndex fineGUIndex(grid_->levelGridView(i+2), gridDim);
           LevelGUIndex coarseGUIndex(grid_->levelGridView(i+1), gridDim);
 
+          typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LevelGridView, Dune::MCMGVertexLayout> LevelLocalMapper;
+          LevelLocalMapper fineLevelLocalMapper(grid_->levelGridView(i+2));
+          LevelLocalMapper coarseLevelLocalMapper(grid_->levelGridView(i+1));
+
           typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
           MatrixCommunicator<LevelGUIndex,
-                             typename GridType::LevelGridView::CollectiveCommunication,
-                             TransferOperatorType> matrixComm(fineGUIndex, coarseGUIndex, grid_->levelGridView(i+1).comm(), 0);
+                             typename GridType::LevelGridView,
+                             TransferOperatorType,
+                             LevelLocalMapper,
+                             LevelLocalMapper> matrixComm(fineGUIndex, coarseGUIndex, grid_->levelGridView(i+1), fineLevelLocalMapper, coarseLevelLocalMapper, 0);
 
           mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>;
           Dune::shared_ptr<TransferOperatorType> transferOperatorMatrix = Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(newTransferOp->getMatrix()));
@@ -302,9 +321,17 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
     VectorCommunicator<GUIndex, typename GridType::LeafGridView::CollectiveCommunication, CorrectionType> vectorComm(*guIndex_,
                                                                                                                      grid_->leafGridView().comm(),
                                                                                                                      0);
-    MatrixCommunicator<GUIndex, typename GridType::LeafGridView::CollectiveCommunication, MatrixType> matrixComm(*guIndex_,
-                                                                                                                 grid_->leafGridView().comm(),
-                                                                                                                 0);
+
+    LocalMapper localMapper(grid_->leafGridView());
+    MatrixCommunicator<GUIndex,
+                       typename GridType::LeafGridView,
+                       MatrixType,
+                       LocalMapper,
+                       LocalMapper> matrixComm(*guIndex_,
+                                               grid_->leafGridView(),
+                                               localMapper,
+                                               localMapper,
+                                               0);
 
     for (int i=0; i<maxTrustRegionSteps_; i++) {
 

dune/gfe/riemanniantrsolver.hh

@@ -51,8 +51,10 @@ class RiemannianTrustRegionSolver
     typedef P3NodalBasis<typename GridType::LeafGridView,double> BasisType;
 #elif defined SECOND_ORDER
     typedef P2NodalBasis<typename GridType::LeafGridView,double> BasisType;
+    typedef P2BasisMapper<typename GridType::LeafGridView> LocalMapper;
 #else
     typedef P1NodalBasis<typename GridType::LeafGridView,double> BasisType;
+    typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView, Dune::VertexLayout> LocalMapper;
 #endif
 
 public: