couplingmanager_binary_cvfe.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_BOUNDARY_COUPLINGMANAGER_BINARY_CVFE_HH
#define DUMUX_MULTIDOMAIN_BOUNDARY_COUPLINGMANAGER_BINARY_CVFE_HH
 
#include <utility>
#include <memory>
#include <tuple>
#include <vector>
#include <unordered_map>
#include <cassert>
 
#include <dune/common/exceptions.hh>
#include <dumux/common/concepts/localdofs_.hh>
#include <dumux/common/concepts/variables_.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/typetraits/griddiscretization.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/cvfe/interpolationpointdata.hh>
#include <dumux/common/typetraits/localdofs_.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dumux/multidomain/couplingmanager.hh>
 
#include <dumux/multidomain/boundary/couplingmapper_boundaryfaces.hh>
 
namespace Dumux {

template<class MDTraits>
class CouplingManagerBinaryCvfe
: public CouplingManager<MDTraits>
{
    using ParentType = CouplingManager<MDTraits>;
 
public:
    using SolutionVectorStorage = typename ParentType::SolutionVectorStorage;
private:
    static_assert(MDTraits::numSubDomains == 2, "CouplingManagerBinaryCvfe requires exactly two subdomains");
 
    template<std::size_t i>
    static constexpr auto domainIdx()
    { return typename MDTraits::template SubDomain<i>::Index{}; }
 
    // the sub domain type tags
    template<std::size_t id>
    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
    template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView;
    template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
    template<std::size_t id> using GridDiscretization = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
 
    static_assert(DiscretizationMethods::isCVFE<std::decay_t<decltype(GridDiscretization<0>::discMethod)>>,
        "CouplingManagerBinaryCvfe requires subdomain 0 to be discretized with a CVFE scheme");
    static_assert(DiscretizationMethods::isCVFE<std::decay_t<decltype(GridDiscretization<1>::discMethod)>>,
        "CouplingManagerBinaryCvfe requires subdomain 1 to be discretized with a CVFE scheme");
 
    template<std::size_t id> using ElementDiscretization = typename GridDiscretization<id>::LocalView;
    template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
    template<std::size_t id> using GridVariablesCache = Concept::GridVariablesCache_t<GridVariables<id>>;
    template<std::size_t id> using Variables = Concept::Variables_t<GridVariables<id>>;
    template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
    template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
 
    using GridVariablesTuple = typename MDTraits::template TupleOfSharedPtr<GridVariables>;
 
    using CouplingMapper = MapperCoupledMatchingBoundaryFaces<MDTraits, CouplingManagerBinaryCvfe<MDTraits>>;
    using CouplingStencil = typename CouplingMapper::CouplingStencil;
 
    template<class ElementSolution, std::size_t domainIdx, std::size_t coupledDomainIdx>
    struct CouplingContext
    {
        ElementDiscretization<coupledDomainIdx> elemDisc;
        const Problem<coupledDomainIdx>& problem;
        ElementSolution elemSol;
 
        template<class ScvOrLocalDof>
        auto vars(const ScvOrLocalDof& scvOrLocalDof) const
        {
            Variables<coupledDomainIdx> variables;
            if constexpr (Concept::FVGridVariables<GridVariables<coupledDomainIdx>>)
            {
                // the variables are still defined per sub-control volume
                if constexpr (Concept::LocalDof<ScvOrLocalDof>)
                    variables.update(elemSol, problem, elemDisc.element(), elemDisc.scv(scvOrLocalDof.index()));
                else
                    variables.update(elemSol, problem, elemDisc.element(), scvOrLocalDof);
            }
            else
                variables.update(elemSol, problem, elemDisc, ipData(elemDisc, scvOrLocalDof));
            return variables;
        }
 
        template<class ScvOrLocalDof>
        auto operator[](const ScvOrLocalDof& scvOrLocalDof) const
        {
            return vars(scvOrLocalDof);
        }
    };
 
public:

    // \{

    void init(std::shared_ptr<Problem<domainIdx<0>()>> domain0Problem,
              std::shared_ptr<Problem<domainIdx<1>()>> domain1Problem,
              GridVariablesTuple&& gridVariables,
              const SolutionVectorStorage& curSol)
    {
        this->setSubProblems(std::make_tuple(domain0Problem, domain1Problem));
        gridVariables_ = std::move(gridVariables);
        this->attachSolution(curSol);
        couplingMapper_.update(*this);
    }
 
    // \}
 

    // \{

    template<std::size_t i, class Assembler>
    void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler) const
    {}

    template<std::size_t i, std::size_t j, class LocalAssemblerI>
    void updateCouplingContext(Dune::index_constant<i> domainI,
                               const LocalAssemblerI& localAssemblerI,
                               Dune::index_constant<j> domainJ,
                               std::size_t dofIdxGlobalJ,
                               const PrimaryVariables<j>& priVarsJ,
                               int pvIdxJ)
    {
        this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
        // couplingContext() recomputes variables on-the-fly via makeCouplingContext_
        // using the correct element from the coupling mapper, so no cache update needed here.
    }
 
    // \}

    template<std::size_t i, class FaceIpData>
    auto couplingContext(Dune::index_constant<i> domainI,
                         const ElementDiscretization<i>& elemDisc,
                         const FaceIpData& faceIpData) const
    { return makeCouplingContext_(domainI, elemDisc, faceIpData); }

    // \{

    template<std::size_t i, std::size_t j>
    const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI,
                                           const Element<i>& element,
                                           Dune::index_constant<j> domainJ) const
    {
        const auto eIdx = Dumux::gridDiscretization(this->problem(domainI)).elementMapper().index(element);
        return couplingMapper_.couplingStencil(domainI, eIdx, domainJ);
    }
 
    // \}

    template<std::size_t i>
    bool isCoupled(Dune::index_constant<i> domainI,
                   const ElementDiscretization<i>& elemDisc,
                   const typename ElementDiscretization<i>::SubControlVolumeFace& scvf) const
    {
        return isCoupled(domainI, elemDisc, elemDisc.intersectionIndex(scvf));
    }

    template<std::size_t i>
    bool isCoupled(Dune::index_constant<i> domainI,
                   const typename ElementDiscretization<i>::SubControlVolume& scv) const
    { return couplingMapper_.isCoupledDof(domainI, scv.dofIndex()); }

    template<std::size_t i>
    bool isCoupled(Dune::index_constant<i> domainI,
                   const ElementDiscretization<i>& elemDisc,
                   std::size_t iIdx) const
    {
        const auto eIdx = Dumux::gridDiscretization(this->problem(domainI)).elementMapper().index(elemDisc.element());
        return couplingMapper_.isCoupled(domainI, eIdx, iIdx);
    }

    template<std::size_t i>
    bool isCoupled(Dune::index_constant<i> domainI,
                   std::size_t dofIdx) const
    { return couplingMapper_.isCoupledDof(domainI, dofIdx); }
 
private:
    template<class ElementDiscretization, class FaceIpData>
    std::size_t intersectionIndex_(const ElementDiscretization& elemDisc,
                                   const FaceIpData& faceIpData) const
    {
        if constexpr (requires { faceIpData.scvfIndex(); })
        {
            const auto& scvf = elemDisc.scvf(faceIpData.scvfIndex());
            return elemDisc.intersectionIndex(scvf);
        }
        else if constexpr (requires { faceIpData.boundaryFaceIndex(); })
            return elemDisc.boundaryFace(faceIpData.boundaryFaceIndex()).intersectionIndex();
        else
            DUNE_THROW(Dune::InvalidStateException, "FaceIpData must provide either scvfIndex() or boundaryFaceIndex()");
    }
 
    template<std::size_t i, class FaceIpData>
    auto makeCouplingContext_(Dune::index_constant<i> domainI,
                              const ElementDiscretization<i>& elemDisc,
                              const FaceIpData& faceIpData) const
    {
        const auto eIdx = Dumux::gridDiscretization(this->problem(domainI)).elementMapper().index(elemDisc.element());
        const auto iIdx = intersectionIndex_(elemDisc, faceIpData);
 
        if (!couplingMapper_.isCoupled(domainI, eIdx, iIdx))
            DUNE_THROW(Dune::InvalidStateException, "No coupling context found at element " << eIdx << " with intersection index " << iIdx);
 
        const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, eIdx, iIdx);
        constexpr auto domainJ = Dune::index_constant<1-i>();
        const auto& otherGridDiscretization = Dumux::gridDiscretization(this->problem(domainJ));
        const auto& otherElement = otherGridDiscretization.element(otherElementIdx);
        auto otherElemDisc = localView(otherGridDiscretization).bind(otherElement);
        const auto elemSol = elementSolution(otherElement, this->curSol(domainJ), otherGridDiscretization);
        return CouplingContext<decltype(elemSol), i, 1-i>{ std::move(otherElemDisc), this->problem(domainJ), std::move(elemSol) };
    }
 
    CouplingMapper couplingMapper_;
    GridVariablesTuple gridVariables_;
};
 
} // end namespace Dumux
 
#endif