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:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
 
#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 ()
    {
        // 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;
    }
 
private:
    // GridView is needed to keep the IndexSet valid
    const 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