diff --git a/dune/gfe/parallel/matrixcommunicator.hh b/dune/gfe/parallel/matrixcommunicator.hh index c934286dcd1eb828b9148fb52281233f0a3bee65..3ddd3b15cf6065616aa0461423bf7facbfe4c0fc 100644 --- a/dune/gfe/parallel/matrixcommunicator.hh +++ b/dune/gfe/parallel/matrixcommunicator.hh @@ -9,7 +9,7 @@ #include <dune/gfe/parallel/mpifunctions.hh> -template<typename RowGlobalMapper, typename GridView, typename MatrixType, typename LocalMapper1, typename LocalMapper2, typename ColumnGlobalMapper=RowGlobalMapper> +template<typename RowGlobalMapper, typename GridView1, typename GridView2, typename MatrixType, typename LocalMapper1, typename LocalMapper2, typename ColumnGlobalMapper=RowGlobalMapper> class MatrixCommunicator { struct TransferMatrixTuple { @@ -45,7 +45,7 @@ class MatrixCommunicator { } public: - MatrixCommunicator(const RowGlobalMapper& rowGlobalMapper, const GridView& gridView, const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root) + MatrixCommunicator(const RowGlobalMapper& rowGlobalMapper, const GridView1& gridView, const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root) : rowGlobalMapper_(rowGlobalMapper), columnGlobalMapper_(rowGlobalMapper), localMapper1_(localMapper1), @@ -53,18 +53,22 @@ public: communicator_(gridView.comm()), root_rank(root) { - setLocalToGlobal(gridView); + setLocalToGlobalRows(gridView); + setLocalToGlobalColumns(gridView); } - MatrixCommunicator(const RowGlobalMapper& rowGlobalMapper, const ColumnGlobalMapper& columnGlobalMapper, const GridView& gridView, const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root) + MatrixCommunicator(const RowGlobalMapper& rowGlobalMapper, const ColumnGlobalMapper& columnGlobalMapper, + const GridView1& gridView1, const GridView2& gridView2, + const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root) : rowGlobalMapper_(rowGlobalMapper), columnGlobalMapper_(columnGlobalMapper), localMapper1_(localMapper1), localMapper2_(localMapper2), - communicator_(gridView.comm()), + communicator_(gridView1.comm()), root_rank(root) { - setLocalToGlobal(gridView); + setLocalToGlobalRows(gridView1); + setLocalToGlobalColumns(gridView2); } MatrixType reduceAdd(const MatrixType& local) @@ -119,27 +123,38 @@ public: private: - void setLocalToGlobal(const GridView& gridView) + void setLocalToGlobalRows(const GridView1& 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 (int codim = 0; codim <= GridView1::dimension; codim++) for (size_t i=0; i<it->subEntities(codim); i++) { - typename RowGlobalMapper::Index localIdx; + typename LocalMapper1::Index localIdx; typename RowGlobalMapper::Index globalIdx; if (localMapper1_.contains(*it,i,codim,localIdx) and rowGlobalMapper_.contains(*it,i,codim,globalIdx)) localToGlobal1_[localIdx] = globalIdx; + } + + } + + void setLocalToGlobalColumns(const GridView2& gridView) + { + localToGlobal2_.resize(localMapper2_.size()); + for (auto it = gridView.template begin<0>(); it != gridView.template end<0>(); ++it) + for (int codim = 0; codim <= GridView2::dimension; codim++) + for (size_t i=0; i<it->subEntities(codim); i++) + { + typename LocalMapper2::Index localIdx; + typename ColumnGlobalMapper::Index globalIdx; if (localMapper2_.contains(*it,i,codim,localIdx) and columnGlobalMapper_.contains(*it,i,codim,globalIdx)) localToGlobal2_[localIdx] = globalIdx; } - } // Mappers for the global numbering @@ -150,7 +165,7 @@ private: const LocalMapper1& localMapper1_; const LocalMapper2& localMapper2_; - const typename GridView::CollectiveCommunication& communicator_; + const typename GridView1::CollectiveCommunication& communicator_; int root_rank; std::vector<typename RowGlobalMapper::Index> localToGlobal1_; diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc index 4a3f567dc8cf3396eb2558b4c4405a7dcb7580c2..8f9d44b4811dbbc8fbefcc1800fb99fb5e5066c8 100644 --- a/dune/gfe/riemanniantrsolver.cc +++ b/dune/gfe/riemanniantrsolver.cc @@ -131,6 +131,7 @@ setup(const GridType& grid, LocalMapper localMapper(grid_->leafGridView()); MatrixCommunicator<GlobalMapper, + typename GridType::LeafGridView, typename GridType::LeafGridView, ScalarMatrixType, LocalMapper, @@ -199,11 +200,12 @@ setup(const GridType& grid, typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType; MatrixCommunicator<GlobalMapper, + typename GridType::LeafGridView, typename GridType::LeafGridView, TransferOperatorType, LocalMapper, LeafP1LocalMapper, - GlobalLeafP1Mapper> matrixComm(*globalMapper_, p1Index, grid_->leafGridView(), localMapper, leafP1LocalMapper, 0); + GlobalLeafP1Mapper> matrixComm(*globalMapper_, p1Index, grid_->leafGridView(), grid_->leafGridView(), localMapper, leafP1LocalMapper, 0); mmgStep->mgTransfer_.back() = new PKtoP1MGTransfer<CorrectionType>; Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(topTransferOp->getMatrix())); @@ -226,10 +228,11 @@ setup(const GridType& grid, typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType; MatrixCommunicator<GlobalLevelP1Mapper, + typename GridType::LevelGridView, typename GridType::LevelGridView, TransferOperatorType, LevelLocalMapper, - LevelLocalMapper> matrixComm(fineGUIndex, coarseGUIndex, grid_->levelGridView(i+1), fineLevelLocalMapper, coarseLevelLocalMapper, 0); + LevelLocalMapper> matrixComm(fineGUIndex, coarseGUIndex, grid_->levelGridView(i+2), 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())); @@ -337,6 +340,7 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve() LocalMapper localMapper(grid_->leafGridView()); MatrixCommunicator<GlobalMapper, + typename GridType::LeafGridView, typename GridType::LeafGridView, MatrixType, LocalMapper,