diff --git a/dune/gfe/parallel/globalindex.hh b/dune/gfe/parallel/globalindex.hh index 80d11aa1520584fdde67c5bdb537fa5a7ef25798..d3c665ed416d5df62360095349faae7f23ee1539 100644 --- a/dune/gfe/parallel/globalindex.hh +++ b/dune/gfe/parallel/globalindex.hh @@ -340,4 +340,108 @@ protected: IndexMap globalLocalMap_; }; + +namespace Dune { + + template <class GridView> + class GlobalP2Mapper + { + public: + + typedef std::map<int,int> IndexMap; + + GlobalP2Mapper(const GridView& gridView) + : gridView_(gridView) + { + static_assert(GridView::dimension==2, "Only implemented for two-dimensional grids"); + + P2BasisMapper<GridView> p2Mapper(gridView); + + GlobalUniqueIndex<GridView,2> globalVertexIndex(gridView); + GlobalUniqueIndex<GridView,1> globalEdgeIndex(gridView); + GlobalUniqueIndex<GridView,0> globalElementIndex(gridView); + + // total number of degrees of freedom + nGlobalEntity_ = globalVertexIndex.nGlobalEntity() + globalEdgeIndex.nGlobalEntity() + globalElementIndex.nGlobalEntity(); + nOwnedLocalEntity_ = globalVertexIndex.nOwnedLocalEntity() + globalEdgeIndex.nOwnedLocalEntity() + globalElementIndex.nOwnedLocalEntity(); + + // Determine + for (auto it = gridView.template begin<0>(); it != gridView.template end<0>(); ++it) + { + // Loop over all vertices + for (size_t i=0; i<it->template count<2>(); i++) + { + //int localIndex = globalVertexIndex.localIndex (*it->template subEntity<2>(i)); + int localIndex = p2Mapper.map(*it, i, 2); + int globalIndex = globalVertexIndex.globalIndex(*it->template subEntity<2>(i)); + + localGlobalMap_[localIndex] = globalIndex; + globalLocalMap_[globalIndex] = localIndex; + } + + // Loop over all edges + for (size_t i=0; i<it->template count<1>(); i++) + { + //int localIndex = globalEdgeIndex.localIndex (*it->template subEntity<1>(i)) + gridView.size(2); + int localIndex = p2Mapper.map(*it, i, 1); + int globalIndex = globalEdgeIndex.globalIndex(*it->template subEntity<1>(i)) + globalVertexIndex.nGlobalEntity(); + + localGlobalMap_[localIndex] = globalIndex; + globalLocalMap_[globalIndex] = localIndex; + } + + // One element degree of freedom for quadrilaterals + if (not it->type().isQuadrilateral()) + DUNE_THROW(Dune::NotImplemented, "for non-quadrilaterals"); + + if (it->type().isQuadrilateral()) + { + //int localIndex = globalEdgeIndex.localIndex (*it->template subEntity<1>(i)) + gridView.size(2); + int localIndex = p2Mapper.map(*it, 0, 0); + int globalIndex = globalElementIndex.globalIndex(*it->template subEntity<0>(0)) + + globalEdgeIndex.nGlobalEntity() + + globalVertexIndex.nGlobalEntity(); + + localGlobalMap_[localIndex] = globalIndex; + globalLocalMap_[globalIndex] = localIndex; + } + + } + + } + + /** \brief Given a local index, retrieve its index globally unique over all processes. */ + int globalIndex(const int& localIndex) const { + return localGlobalMap_.find(localIndex)->second; + } + + int localIndex(const int& globalIndex) const { + return globalLocalMap_.find(globalIndex)->second; + } + + unsigned int nGlobalEntity() const + { + return nGlobalEntity_; + } + + unsigned int nOwnedLocalEntity() const + { + return nOwnedLocalEntity_; + } + + const GridView& getGridView() const { + return gridView_; + } + + const GridView gridView_; + + IndexMap localGlobalMap_; + IndexMap globalLocalMap_; + + size_t nOwnedLocalEntity_; + size_t nGlobalEntity_; + + }; + +} #endif /* GLOBALUNIQUEINDEX_HH_ */