couplingmapperbase.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_MAPPER_BASE_HH
#define DUMUX_FACETCOUPLING_MAPPER_BASE_HH
 
#include <vector>
#include <unordered_map>
#include <algorithm>
#include <memory>
 
#include <dune/common/indices.hh>
#include <dune/common/exceptions.hh>
#include <dumux/common/indextraits.hh>
 
namespace Dumux {
 
template< class BulkFVG,
          class LowDimFVG,
          std::size_t bulkId,
          std::size_t lowDimId>
class FacetCouplingMapperBase
{
    using BulkGridView = typename BulkFVG::GridView;
    using BulkIndexType = typename IndexTraits<BulkGridView>::GridIndex;
 
    using LowDimGridView = typename LowDimFVG::GridView;
    using LowDimElement = typename LowDimGridView::template Codim<0>::Entity;
    using LowDimIndexType = typename IndexTraits<LowDimGridView>::GridIndex;
 
    // make sure the grid geometry combination makes sense
    static constexpr int bulkDim = BulkGridView::dimension;
    static constexpr int lowDimDim = LowDimGridView::dimension;
    static_assert(bulkDim == lowDimDim+1, "Lower-dimensional geometry is not of codim 1 w.r.t. bulk geometry!");
 
    // discretization method of the bulk domain
    static constexpr auto bulkDiscMethod = BulkFVG::discMethod;
 
    // helper struct to check validity of the given domain id offset
    template< class Embeddings >
    class IsValidDomainId
    {
        using GCBulkGridView = typename Embeddings::template GridView<bulkId>;
        using GCLowDimGridView = typename Embeddings::template GridView<lowDimId>;
        static constexpr bool bulkMatch = std::is_same<GCBulkGridView, BulkGridView>::value;
        static constexpr bool lowDimMatch = std::is_same<GCLowDimGridView, LowDimGridView>::value;
        static_assert(bulkMatch, "The bulk domain id does not match the provided bulk grid geometry");
        static_assert(lowDimMatch, "The lowDim domain id does not match the provided lowDim grid geometry");
 
    public:
        static constexpr bool value = bulkMatch && lowDimMatch;
    };
 
    struct BulkCouplingData
    {
        std::vector< LowDimIndexType > couplingStencil;
 
        std::vector< LowDimIndexType > couplingElementStencil;
 
        std::unordered_map< LowDimIndexType, std::vector<BulkIndexType> > dofToCouplingScvfMap;
 
        std::unordered_map< LowDimIndexType, std::vector<BulkIndexType> > elementToScvfMap;
    };
 
    struct LowDimCouplingData
    {
        std::vector< BulkIndexType > couplingStencil;
 
        using Embedment = std::pair< BulkIndexType, std::vector<BulkIndexType> >;
        std::vector< Embedment > embedments;
    };
 
    using BulkCouplingMap = std::unordered_map<BulkIndexType, BulkCouplingData>;
    using LowDimCouplingMap = std::unordered_map<LowDimIndexType, LowDimCouplingData>;
 
    using LowDimStencil = std::vector<LowDimIndexType>;
    using BulkStencil = std::vector<BulkIndexType>;
 
    template<std::size_t id>
    using GridIdType = Dune::index_constant<id>;
 
public:
    static constexpr auto bulkGridId = GridIdType<bulkId>();
    static constexpr auto facetGridId = GridIdType<lowDimId>();
 
    template<std::size_t id>
    using Stencil = typename std::conditional<id == bulkId, BulkStencil, LowDimStencil>::type;
 
    template<std::size_t i, std::size_t j>
    using CouplingMap = typename std::conditional<i == bulkId, BulkCouplingMap, LowDimCouplingMap>::type;
 
    template<int dim>
    static constexpr GridIdType< (dim == bulkDim ? bulkId : lowDimId) > gridId()
    { return GridIdType< (dim == bulkDim ? bulkId : lowDimId) >(); }
 
    template< class Embeddings >
    void update(const BulkFVG& bulkFvGridGeometry,
                const LowDimFVG& lowDimFvGridGeometry,
                std::shared_ptr<const Embeddings> embeddings)
    { DUNE_THROW(Dune::NotImplemented, "Implementation does not provide an update() function."); }
 
    const BulkCouplingMap& couplingMap(GridIdType<bulkId>, GridIdType<lowDimId>) const
    { return bulkCouplingData_; }
 
    const LowDimCouplingMap& couplingMap(GridIdType<lowDimId>, GridIdType<bulkId>) const
    { return lowDimCouplingData_; }
 
protected:
 
    template< class Embeddings, typename AddCouplingEntryPolicy >
    void update_(const BulkFVG& bulkFvGridGeometry,
                 const LowDimFVG& lowDimFvGridGeometry,
                 std::shared_ptr<const Embeddings> embeddings,
                 AddCouplingEntryPolicy&& addCouplingEntryPolicy)
    {
        // some static assertions on the grid creator
        static_assert(IsValidDomainId<Embeddings>::value, "Grid type mismatch. Please review the provided domain id offset.");
 
        // clear data
        bulkCouplingData_.clear();
        lowDimCouplingData_.clear();
 
        // set up maps between element indices and insertion indices
        const auto bulkInsertionToElemIdxMap = makeInsertionToGridIndexMap_(embeddings, bulkFvGridGeometry);
 
        // set up coupling maps coming from the low dim domain
        for (const auto& element : elements(lowDimFvGridGeometry.gridView()))
        {
            auto adjoinedEntities = embeddings->template adjoinedEntityIndices<lowDimId>(element);
 
            // proceed only if element is embedded
            if (adjoinedEntities.size() == 0)
                continue;
 
            // turn (insertion) indices into actual grid element indices ...
            std::for_each(adjoinedEntities.begin(), adjoinedEntities.end(), [&] (auto& idx) { idx = bulkInsertionToElemIdxMap.at(idx); });
 
            // ... and add them
            addCouplingEntryPolicy(std::move(adjoinedEntities), element, lowDimFvGridGeometry, bulkFvGridGeometry);
        }
    }
 
    template< class GridGeometry>
    std::vector< typename IndexTraits<typename GridGeometry::GridView>::GridIndex >
    extractNodalDofs_(const typename GridGeometry::GridView::template Codim<0>::Entity& element,
                      const GridGeometry& gridGeometry)
    {
        static constexpr int dim = GridGeometry::GridView::dimension;
        using GridIndexType = typename IndexTraits<typename GridGeometry::GridView>::GridIndex;
 
        const auto numCorners = element.subEntities(dim);
        std::vector< GridIndexType > nodalDofs(numCorners);
        for (unsigned int i = 0; i < numCorners; ++i)
            nodalDofs[i] = gridGeometry.vertexMapper().subIndex(element, i, dim);
 
        return nodalDofs;
    }
 
    BulkCouplingMap& couplingMap_(GridIdType<bulkId>, GridIdType<lowDimId>)
    { return bulkCouplingData_; }
 
    LowDimCouplingMap& couplingMap_(GridIdType<lowDimId>, GridIdType<bulkId>)
    { return lowDimCouplingData_; }
 
private:
 
    template< class Embeddings, class GridGeometry>
    std::unordered_map< typename IndexTraits<typename GridGeometry::GridView>::GridIndex, typename IndexTraits<typename GridGeometry::GridView>::GridIndex >
    makeInsertionToGridIndexMap_(std::shared_ptr<const Embeddings> embeddings, const GridGeometry& gridGeometry) const
    {
        using GridIndexType = typename IndexTraits<typename GridGeometry::GridView>::GridIndex;
 
        std::unordered_map< GridIndexType, GridIndexType > map;
        for (const auto& e : elements(gridGeometry.gridView()))
            map.insert( std::make_pair( embeddings->template insertionIndex<bulkId>(e), gridGeometry.elementMapper().index(e) ) );
 
        return map;
    }
 
    BulkCouplingMap bulkCouplingData_;     
    LowDimCouplingMap lowDimCouplingData_; 
};
 
} // end namespace Dumux
 
#endif