reorderingdofmapper.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
 
#ifndef DUMUX_COMMON_REORDERING_DOF_MAPPER_HH
#define DUMUX_COMMON_REORDERING_DOF_MAPPER_HH
 
#if HAVE_PTSCOTCH
 
#include <dune/common/timer.hh>
#include <dune/grid/common/mapper.hh>
 
#include <dumux/linear/scotchbackend.hh>
 
namespace Dumux {
 
template<class GridView>
class ReorderingDofMapper
: public Dune::Mapper<typename GridView::Grid, ReorderingDofMapper<GridView>, typename GridView::IndexSet::IndexType>
{
    using Index = typename GridView::IndexSet::IndexType;
    using ParentType = typename Dune::Mapper<typename GridView::Grid, ReorderingDofMapper<GridView>, Index>;
    using Element = typename GridView::template Codim<0>::Entity;
public:
 
    template<class Layout>
    ReorderingDofMapper (const GridView& gridView, Layout&& layout)
    : gridView_(gridView)
    , indexSet_(&gridView.indexSet())
    , codimension_(layout(indexSet_->types(0)[0], GridView::dimension) ? 0 : GridView::dimension)
    {
        update_();
    }
 
    template<class EntityType>
    Index index (const EntityType& e) const
    {
        // map the index using the permutation obtained from the reordering algorithm
        return static_cast<Index>(permutation_[indexSet_->index(e)]);
    }
 
    Index subIndex (const Element& e, int i, unsigned int codim) const
    {
        return indexSet_->subIndex(e, i, codim);
    }
 
    std::size_t size () const
    {
        return indexSet_->size(codimension_);
    }
 
    template<class EntityType>
    bool contains (const EntityType& e, Index& result) const
    {
        result = index(e);
        return true;
    }
 
    bool contains (const Element& e, int i, int cc, Index& result) const
    {
        result = indexSet_->subIndex(e, i, cc);
        return true;
    }
 
    void update (const GridView& gridView)
    {
        gridView_ = gridView;
        indexSet_ = &gridView_.indexSet();
        update_();
    }
 
    void update (GridView&& gridView)
    {
        gridView_ = std::move(gridView);
        indexSet_ = &gridView_.indexSet();
        update_();
    }
 
private:
    void update_()
    {
        // Compute scotch reordering
        Dune::Timer watch;
        // Create the graph as an adjacency list
        std::vector<std::vector<Index>> graph(size());
 
        // dofs on element centers (cell-centered methods)
        if (codimension_ == 0)
        {
            for (const auto& element : elements(gridView_))
            {
                auto eIdx = indexSet_->index(element);
                for (const auto& intersection : intersections(gridView_, element))
                {
                    if (intersection.neighbor())
                        graph[eIdx].push_back(indexSet_->index(intersection.outside()));
                }
            }
        }
        // dof on vertices (box method)
        else
        {
            for (const auto& element : elements(gridView_))
            {
                auto eIdx = indexSet_->index(element);
                for (int vIdxLocal = 0; vIdxLocal < element.subEntities(codimension_); ++vIdxLocal)
                {
                    auto vIdxGlobal = indexSet_->subIndex(element, vIdxLocal, codimension_);
                    graph[vIdxGlobal].push_back(eIdx);
                }
            }
        }
 
        permutation_ = ScotchBackend<Index>::computeGPSReordering(graph);
        std::cout << "Scotch backend reordered index set of size " << size()
                  << " in " << watch.elapsed() << " seconds." << std::endl;
    }
    // GridView is needed to keep the IndexSet valid
    GridView gridView_;
    const typename GridView::IndexSet* indexSet_;
    const int codimension_;
    // the map resulting from the reordering
    std::vector<int> permutation_;
};
 
} // end namespace Dumux
 
#else
 
#warning "PTSCOTCH was not found on your system. Dumux::ReorderingDofMapper needs it to work -> fallback to MCMGMapper without reordering!"
#include <dune/grid/common/mcmgmapper.hh>
namespace Dumux {
template<class GridView>
using ReorderingDofMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
} // end namespace Dumux
 
#endif // HAVE_PTSCOTCH
#endif