multidomain/facet/gridmanager.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_FACETCOUPLING_GRID_MANAGER_HH
#define DUMUX_FACETCOUPLING_GRID_MANAGER_HH
 
#include <cassert>
#include <map>
#include <tuple>
#include <type_traits>
#include <vector>
 
#include <dune/common/exceptions.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/grid/common/gridfactory.hh>
 
#include <dumux/io/grid/griddata.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/indextraits.hh>
#include <dumux/common/typetraits/utility.hh>
 
#include "gmshreader.hh"
 
namespace Dumux {
namespace FCGridManagerChecks {
 
    // The grid creator and grid data classes provided below
    // require that the grids passed as template arguments are
    // ordered in descending grid dimension and they must have
    // the same world dimension. The following helper structs
    // verify these conditions and are reused below.
    // evaluates if dim is descending or dimworld is equal for two grids
    template<bool checkDimWorld, typename G1, typename G2>
    static constexpr bool evalCondition()
    {
        return checkDimWorld ? int(G1::dimensionworld) == int(G2::dimensionworld)
                             : int(G1::dimension) > int(G2::dimension);
    }
 
    // helper structs to evaluate conditions on the grid
    template<bool checkDimWorld, typename... Gs> struct FulfillConditions;
    template<bool checkDimWorld, typename G1, typename... Gs>
    struct FulfillConditions<checkDimWorld, G1, Gs...>
    {
        using G2 = typename std::tuple_element_t<0, std::tuple<Gs...>>;
        static constexpr bool value = evalCondition<checkDimWorld, G1, G2>() && FulfillConditions<checkDimWorld, Gs...>::value;
    };
    template<bool checkDimWorld, typename G1, typename G2>
    struct FulfillConditions<checkDimWorld, G1, G2> { static constexpr bool value = evalCondition<checkDimWorld, G1, G2>(); };
}
 
template<typename... Grids>
class FacetCouplingGridDataWrapper
{
    // make sure all grids have the same world dimension and are ordered in descending dimension
    static_assert(FCGridManagerChecks::FulfillConditions<false, Grids...>::value, "All grids must have the same world dimension!");
    static_assert(FCGridManagerChecks::FulfillConditions<true, Grids...>::value, "Grids must be ordered w.r.t the dimension in descending order!");
 
    static constexpr std::size_t numGrids = sizeof...(Grids);
    template<std::size_t id> using Grid = typename std::tuple_element_t<id, std::tuple<Grids...>>;
    template<std::size_t id> using GridDataPtr = std::shared_ptr< GridData<Grid<id>> >;
    template<std::size_t id> using Intersection = typename Grid<id>::LeafGridView::Intersection;
 
public:
    template<std::size_t id> using GridData = GridData<Grid<id>>;
 
    template<std::size_t id>
    void setGridData(GridData<id>&& gridData)
    {
        std::get<id>(gridDataPtrTuple_) = std::make_shared<GridData<id>>( std::move(gridData) );
    }
 
    template<std::size_t id>
    std::shared_ptr<const GridData<id>> getSubDomainGridData() const
    { return std::get<id>(gridDataPtrTuple_); }
 
    template<std::size_t id>
    int getElementDomainMarker(const typename Grid<id>::template Codim<0>::Entity& element) const
    { return std::get<id>(gridDataPtrTuple_)->getElementDomainMarker(element); }
 
    template<std::size_t id>
    int getBoundaryDomainMarker(const typename Grid<id>::LeafGridView::Intersection& is) const
    { return std::get<id>(gridDataPtrTuple_)->getBoundaryDomainMarker(is); }
 
    template<std::size_t id>
    int getBoundaryDomainMarker(int boundarySegmentIndex) const
    { return std::get<id>(gridDataPtrTuple_)->getBoundaryDomainMarker(boundarySegmentIndex); }
 
    template<std::size_t id>
    bool wasInserted(const Intersection<id>& intersection) const
    { return std::get<id>(gridDataPtrTuple_)->wasInserted(intersection); }
 
private:
    using GridDataPtrTuple = typename makeFromIndexedType<std::tuple, GridDataPtr, std::make_index_sequence<numGrids>>::type;
    GridDataPtrTuple gridDataPtrTuple_;
};
 
template<typename... Grids>
class FacetCouplingEmbeddings
{
    // make sure all grids have the same world dimension and are ordered in descending dimension
    static_assert(FCGridManagerChecks::FulfillConditions<false, Grids...>::value, "All grids must have the same world dimension!");
    static_assert(FCGridManagerChecks::FulfillConditions<true, Grids...>::value, "Grids must be ordered w.r.t the dimension in descending order!");
 
    template<std::size_t id> using Grid = typename std::tuple_element_t<id, std::tuple<Grids...>>;
    template<std::size_t id> using GridFactory = typename Dune::GridFactory<Grid<id>>;
 
    using GIType = typename IndexTraits< typename Grid<0>::LeafGridView >::GridIndex;
    using EmbedmentMap = std::unordered_map<GIType, std::vector<GIType>>;
 
public:
    template<std::size_t id> using GridView = typename Grid<id>::LeafGridView;
 
    static constexpr std::size_t numGrids = sizeof...(Grids);
    static constexpr int bulkGridId = 0;
    static constexpr int bulkDim = Grid<bulkGridId>::dimension;
 
    using BulkGridView = GridView<bulkGridId>;
    using GridIndexType = GIType;
 
    template<std::size_t id>
    const GridView<id>& gridView() const
    { return *std::get<id>(gridViewPtrTuple_); }
 
    template<std::size_t id, class Entity>
    GridIndexType insertionIndex(const Entity& entity) const
    { return std::get<id>(gridFactoryPtrTuple_)->insertionIndex(entity); }
 
    template<std::size_t id>
    typename std::unordered_map< GridIndexType, std::vector<GridIndexType> >::mapped_type
    embeddedEntityIndices(const typename Grid<id>::template Codim<0>::Entity& element) const
    {
        const auto& map = embeddedEntityMaps_[id];
        auto it = map.find( std::get<id>(gridFactoryPtrTuple_)->insertionIndex(element) );
        if (it != map.end()) return it->second;
        else return typename std::unordered_map< GridIndexType, std::vector<GridIndexType> >::mapped_type();
    }
 
    template<std::size_t id>
    typename std::unordered_map< GridIndexType, std::vector<GridIndexType> >::mapped_type
    adjoinedEntityIndices(const typename Grid<id>::template Codim<0>::Entity& element) const
    {
        const auto& map = adjoinedEntityMaps_[id];
        auto it = map.find( std::get<id>(gridFactoryPtrTuple_)->insertionIndex(element) );
        if (it != map.end()) return it->second;
        else return typename std::unordered_map< GridIndexType, std::vector<GridIndexType> >::mapped_type();
    }
 
    const std::unordered_map< GridIndexType, std::vector<GridIndexType> >& embeddedEntityMap(std::size_t id) const
    { assert(id < numGrids); return embeddedEntityMaps_[id]; }
 
    std::unordered_map< GridIndexType, std::vector<GridIndexType> >& embeddedEntityMap(std::size_t id)
    { assert(id < numGrids); return embeddedEntityMaps_[id]; }
 
    const std::unordered_map< GridIndexType, std::vector<GridIndexType> >& adjoinedEntityMap(std::size_t id) const
    { assert(id < numGrids); return adjoinedEntityMaps_[id]; }
 
    std::unordered_map< GridIndexType, std::vector<GridIndexType> >& adjoinedEntityMap(std::size_t id)
    { assert(id < numGrids); return adjoinedEntityMaps_[id]; }
 
    const std::vector<GridIndexType>& gridHierarchyIndices(std::size_t id) const
    { assert(id < numGrids); return gridVertexIndices_[id]; }
 
    std::size_t numVerticesInHierarchy() const
    { return numVerticesInHierarchy_; }
 
    template<std::size_t id>
    void setData( std::shared_ptr<Grid<id>> gridPtr,
                  std::shared_ptr<GridFactory<id>> gridFactoryPtr,
                  EmbedmentMap&& embeddedEntityMap,
                  EmbedmentMap&& adjoinedEntityMap,
                  std::vector<GridIndexType>&& gridVertexIndices,
                  std::size_t numVerticesInHierarchy )
    {
        std::get<id>(gridViewPtrTuple_) = std::make_shared<GridView<id>>(gridPtr->leafGridView());
        std::get<id>(gridFactoryPtrTuple_) = gridFactoryPtr;
        embeddedEntityMaps_[id] = std::move(embeddedEntityMap);
        adjoinedEntityMaps_[id] = std::move(adjoinedEntityMap);
        gridVertexIndices_[id] = std::move(gridVertexIndices);
        numVerticesInHierarchy_ = numVerticesInHierarchy;
    }
 
private:
    std::array<EmbedmentMap, numGrids> embeddedEntityMaps_;
    std::array<EmbedmentMap, numGrids> adjoinedEntityMaps_;
 
    std::size_t numVerticesInHierarchy_;
    std::array<std::vector<GridIndexType>, numGrids> gridVertexIndices_;
 
    using Indices = std::make_index_sequence<numGrids>;
    template<std::size_t id> using GridViewPtr = std::shared_ptr<GridView<id>>;
    using GridPtrTuple = typename makeFromIndexedType<std::tuple, GridViewPtr, Indices>::type;
    GridPtrTuple gridViewPtrTuple_;
 
    template<std::size_t id> using GridFactoryPtr = std::shared_ptr< Dune::GridFactory<Grid<id>> >;
    using GridFactoryPtrTuple = typename makeFromIndexedType<std::tuple, GridFactoryPtr, Indices>::type;
    GridFactoryPtrTuple gridFactoryPtrTuple_;
};
 
template<typename... Grids>
class FacetCouplingGridManager
{
    // make sure all grids have the same world dimension and are ordered in descending dimension
    static_assert(FCGridManagerChecks::FulfillConditions<false, Grids...>::value, "All grids must have the same world dimension!");
    static_assert(FCGridManagerChecks::FulfillConditions<true, Grids...>::value, "Grids must be ordered w.r.t the dimension in descending order!");
 
    // we use a wrapper class for the grid data containing the data on all grids
    using GridDataWrapper = FacetCouplingGridDataWrapper<Grids...>;
public:
    template<std::size_t id> using Grid = typename std::tuple_element_t<id, std::tuple<Grids...>>;
    template<std::size_t id> using GridPtr = typename std::shared_ptr< Grid<id> >;
 
    static constexpr std::size_t numGrids = sizeof...(Grids);
    static constexpr int bulkGridId = 0;
 
    using GridData = GridDataWrapper;
    using Embeddings = FacetCouplingEmbeddings<Grids...>;
 
    template<std::size_t id>
    const Grid<id>& grid() const
    { return *std::get<id>(gridPtrTuple_); }
 
    std::shared_ptr<const GridData> getGridData() const
    {
        if (!enableEntityMarkers_)
            DUNE_THROW(Dune::IOError, "No grid data available");
        return gridDataPtr_;
    }
 
    std::shared_ptr<const Embeddings> getEmbeddings() const
    { return embeddingsPtr_; }
 
    void init(const std::string& paramGroup = "")
    {
        // reset the grid & embedding data
        gridDataPtr_ = std::make_shared<GridData>();
        embeddingsPtr_ = std::make_shared<Embeddings>();
 
        // get filename and determine grid file extension
        const auto fileName = getParamFromGroup<std::string>(paramGroup, "Grid.File");
        const auto ext = getFileExtension(fileName);
 
        // get some parameters
        const bool verbose = getParamFromGroup<bool>(paramGroup, "Grid.Verbosity", false);
        const bool domainMarkers = getParamFromGroup<bool>(paramGroup, "Grid.DomainMarkers", false);
        const bool boundarySegments = getParamFromGroup<bool>(paramGroup, "Grid.BoundarySegments", false);
 
        // forward to the corresponding reader
        if (ext == "msh")
        {
            const auto thresh = getParamFromGroup<std::size_t>(paramGroup, "Grid.GmshPhysicalEntityThreshold", 0);
            FacetCouplingGmshReader<Grid<bulkGridId>, numGrids> gmshReader;
            gmshReader.read(fileName, (boundarySegments ? thresh : 0), verbose);
            passDataFromReader(gmshReader, domainMarkers, boundarySegments);
        }
        else
            DUNE_THROW(Dune::NotImplemented, "Reader for grid files of type ." + ext);
 
        // find out if entity markers are active
        enableEntityMarkers_ = domainMarkers || boundarySegments;
    }
 
    void loadBalance()
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(gridPtrTuple_)), [&](const auto id)
        {
            std::get<id>(this->gridPtrTuple_)->loadBalance();
        });
    }
 
protected:
    template<std::size_t id>
    Grid<id>& grid_()
    { return *std::get<id>(gridPtrTuple_); }
 
    std::shared_ptr<Embeddings> getEmbeddings_()
    { return embeddingsPtr_; }
 
private:
    static std::string getFileExtension(const std::string& fileName)
    {
        const auto pos = fileName.rfind('.', fileName.length());
        if (pos != std::string::npos)
            return(fileName.substr(pos+1, fileName.length() - pos));
        else
            DUNE_THROW(Dune::IOError, "Please provide an extension for your grid file ('"<< fileName << "')!");
    }
 
    template<typename MeshFileReader>
    void passDataFromReader(MeshFileReader& reader, bool domainMarkers, bool boundarySegments)
    {
        const auto& vertices = reader.gridVertices();
 
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(gridPtrTuple_)), [&](const auto id)
        {
            using GridFactory = Dune::GridFactory<Grid<id>>;
            auto factoryPtr = std::make_shared<GridFactory>();
 
            // insert grid vertices
            for (const auto idx : reader.vertexIndices(id))
                factoryPtr->insertVertex(vertices[idx]);
 
            // insert elements
            for (const auto& e : reader.elementData(id))
                factoryPtr->insertElement(e.gt, e.cornerIndices);
 
            // insert boundary segments
            if (boundarySegments)
                for (const auto& segment : reader.boundarySegmentData(id))
                    factoryPtr->insertBoundarySegment(segment);
 
            // make grid
            auto gridPtr = std::shared_ptr<Grid<id>>(factoryPtr->createGrid());
 
            // maybe create and set grid data object
            if (domainMarkers || boundarySegments)
            {
                typename GridDataWrapper::template GridData<id> gridData( gridPtr,
                                                                          factoryPtr,
                                                                          std::move(reader.elementMarkerMap(id)),
                                                                          std::move(reader.boundaryMarkerMap(id)) );
                gridDataPtr_->template setGridData<id>( std::move(gridData) );
            }
 
            // copy the embeddings
            embeddingsPtr_->template setData<id>( gridPtr,
                                                  factoryPtr,
                                                  std::move(reader.embeddedEntityMap(id)),
                                                  std::move(reader.adjoinedEntityMap(id)),
                                                  std::move(reader.vertexIndices(id)),
                                                  vertices.size() );
 
            // set the grid pointer
            std::get<id>(gridPtrTuple_) = gridPtr;
        });
    }
 
    using Indices = std::make_index_sequence<numGrids>;
    using GridPtrTuple = typename makeFromIndexedType<std::tuple, GridPtr, Indices>::type;
    GridPtrTuple gridPtrTuple_;
 
    bool enableEntityMarkers_;
    std::shared_ptr<GridData> gridDataPtr_;
 
    std::shared_ptr<Embeddings> embeddingsPtr_;
};
 
} // end namespace Dumux
 
#endif