couplingmapper_boundaryfaces.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-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
 
#ifndef DUMUX_MULTIDOMAIN_COUPLINGMAPPER_BOUNDARYFACES_HH
#define DUMUX_MULTIDOMAIN_COUPLINGMAPPER_BOUNDARYFACES_HH
 
#include <algorithm>
#include <array>
#include <iostream>
#include <map>
#include <tuple>
#include <utility>
#include <vector>
 
#include <dune/common/timer.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/indices.hh>
 
#include <dumux/common/typetraits/griddiscretization.hh>
#include <dumux/geometry/intersectingentities.hh>
 
namespace Dumux {

template<class MDTraits, class CouplingManager>
class MapperCoupledMatchingBoundaryFaces
{
    template<std::size_t i> using GridDiscretization = typename MDTraits::template SubDomain<i>::GridGeometry;
    template<std::size_t i> using GridView = typename GridDiscretization<i>::GridView;
 
    template<std::size_t i>
    static constexpr auto domainIdx()
    { return typename MDTraits::template SubDomain<i>::Index{}; }
 
    struct FaceInfo
    {
        std::size_t eIdxOutside;
        std::size_t iIdxOutside;
    };
 
    using FlipIntersectionMapType = std::map<std::pair<std::size_t,std::size_t>, FaceInfo>;
    using MapType = std::map<std::size_t, std::vector<std::size_t>>;
 
    static constexpr std::size_t numSD = MDTraits::numSubDomains;
 
public:
    using CouplingStencil = std::vector<std::size_t>;

    void update(const CouplingManager& couplingManager)
    {
        static_assert(numSD == 2, "More than two subdomains not implemented!");
 
        Dune::Timer watch;
        std::cout << "Initializing the coupling map..." << std::endl;
 
        for (std::size_t domIdx = 0; domIdx < numSD; ++domIdx)
            stencils_[domIdx].clear();
 
        std::get<0>(intersectionInfo_).clear();
        std::get<1>(intersectionInfo_).clear();
 
        const auto& problem0 = couplingManager.problem(domainIdx<0>());
        const auto& problem1 = couplingManager.problem(domainIdx<1>());
        const auto& gridDisc0 = Dumux::gridDiscretization(problem0);
        const auto& gridDisc1 = Dumux::gridDiscretization(problem1);
 
        isCoupledDof_[0].assign(gridDisc0.numDofs(), false);
        isCoupledDof_[1].assign(gridDisc1.numDofs(), false);
 
        auto elemDisc0 = localView(gridDisc0);
        auto elemDisc1 = localView(gridDisc1);
 
        for (const auto& element1 : elements(gridDisc1.gridView()))
        {
            elemDisc1.bindElement(element1);
 
            // the outside elements already added to the stencils for this element,
            // which may touch the coupling interface with more than one intersection
            std::vector<std::size_t> coupledOutsideElements;
 
            for (const auto& boundaryFace1 : boundaryFaces(elemDisc1))
            {
                // get elements intersecting with the face center
                // for robustness add epsilon in normal direction
                const auto eps = boundaryFace1.area()*1e-8;
                auto globalPos = boundaryFace1.center(); globalPos.axpy(eps, boundaryFace1.unitOuterNormal());
                const auto indices = intersectingEntities(globalPos, gridDisc0.boundingBoxTree());
 
                // skip if no intersection was found
                if (indices.empty())
                    continue;
 
                // sanity check
                if (indices.size() > 1)
                    DUNE_THROW(Dune::InvalidStateException, "Are you sure your sub-domain grids are conformingly discretized on the common interface?");
 
                // add the pair to the multimap
                const auto eIdx1 = gridDisc1.elementMapper().index(element1);
                const auto eIdx0 = indices[0];
                const auto& element0 = gridDisc0.element(eIdx0);
                elemDisc0.bindElement(element0);
 
                // add the element dofs once per pair of coupled elements
                if (std::find(coupledOutsideElements.begin(), coupledOutsideElements.end(), eIdx0) == coupledOutsideElements.end())
                {
                    coupledOutsideElements.push_back(eIdx0);
 
                    for(auto&& localDof : localDofs(elemDisc1))
                        stencils_[0][eIdx0].push_back(localDof.dofIndex());
 
                    for(auto&& localDof : localDofs(elemDisc0))
                        stencils_[1][eIdx1].push_back(localDof.dofIndex());
                }
 
                for (const auto& boundaryFace0 : boundaryFaces(elemDisc0))
                {
                    const auto dist = (boundaryFace1.center() - boundaryFace0.center()).two_norm();
                    if (dist > eps)
                        continue;
 
                    for(const auto& localDof : localDofs(elemDisc0, boundaryFace0))
                        isCoupledDof_[0][localDof.dofIndex()] = true;
 
                    for(const auto& localDof : localDofs(elemDisc1, boundaryFace1))
                        isCoupledDof_[1][localDof.dofIndex()] = true;
 
                    std::get<0>(intersectionInfo_)[std::make_pair(eIdx0, static_cast<std::size_t>(boundaryFace0.intersectionIndex()))]
                        = FaceInfo{eIdx1, static_cast<std::size_t>(boundaryFace1.intersectionIndex())};
                    std::get<1>(intersectionInfo_)[std::make_pair(eIdx1, static_cast<std::size_t>(boundaryFace1.intersectionIndex()))]
                        = FaceInfo{eIdx0, static_cast<std::size_t>(boundaryFace0.intersectionIndex())};
                }
            }
        }
 
        // dofs are shared between elements, so the same dof may have been added several times.
        // Duplicate entries would make the coupling derivatives be accumulated more than once.
        for (auto& stencil : stencils_)
            for (auto& [eIdx, dofs] : stencil)
            {
                std::sort(dofs.begin(), dofs.end());
                dofs.erase(std::unique(dofs.begin(), dofs.end()), dofs.end());
            }
 
        std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
    }

    template<std::size_t i, std::size_t j>
    const std::vector<std::size_t>& couplingStencil(Dune::index_constant<i> domainI,
                                                    const std::size_t eIdxI,
                                                    Dune::index_constant<j> domainJ) const
    {
        if (isCoupledElement(domainI, eIdxI))
            return stencils_[i].at(eIdxI);
        else
            return emptyStencil_;
    }

    template<std::size_t i>
    bool isCoupledElement(Dune::index_constant<i>, std::size_t eIdx) const
    {
        return static_cast<bool>(stencils_[i].count(eIdx));
    }

    template<std::size_t i>
    bool isCoupledDof(Dune::index_constant<i>, std::size_t dofIdx) const
    {
        return isCoupledDof_[i][dofIdx];
    }

    template<std::size_t i>
    bool isCoupled(Dune::index_constant<i> domainI,
                   std::size_t eIdxI,
                   std::size_t iIdxI) const
    {
        return static_cast<bool>(std::get<i>(intersectionInfo_).count(std::make_pair(eIdxI, iIdxI)));
    }

    template<std::size_t i>
    std::size_t outsideElementIndex(Dune::index_constant<i> domainI,
                                    std::size_t eIdxI,
                                    std::size_t iIdxI) const
    {
        return std::get<i>(intersectionInfo_).at(std::make_pair(eIdxI, iIdxI)).eIdxOutside;
    }

    template<std::size_t i>
    std::size_t flipIntersectionIndex(Dune::index_constant<i> domainI,
                                      std::size_t eIdxI,
                                      std::size_t iIdxI) const
    {
        return std::get<i>(intersectionInfo_).at(std::make_pair(eIdxI, iIdxI)).iIdxOutside;
    }
 
private:
    std::array<MapType, numSD> stencils_;
    std::vector<std::size_t> emptyStencil_;
    std::array<std::vector<bool>, numSD> isCoupledDof_;
    std::tuple<FlipIntersectionMapType, FlipIntersectionMapType> intersectionInfo_;
};
 
} // end namespace Dumux
 
#endif