multidomain/facet/couplingmanager.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_FACETCOUPLING_MANAGER_HH
#define DUMUX_FACETCOUPLING_MANAGER_HH
 
#include <dumux/common/properties.hh>
#include <dumux/common/indextraits.hh>
 
#include <dumux/discretization/method.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dumux/discretization/evalsolution.hh>
 
namespace Dumux {
namespace FacetCoupling{
 
template<class VolumeVariables, class Problem, class SolutionVector, class FVGeometry>
void makeInterpolatedVolVars(VolumeVariables& volVars,
                             const Problem& problem,
                             const SolutionVector& sol,
                             const FVGeometry& fvGeometry,
                             const typename FVGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
                             const typename FVGeometry::GridGeometry::GridView::template Codim<0>::Entity::Geometry& elemGeom,
                             const typename FVGeometry::GridGeometry::GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate& pos)
{
    // interpolate solution and set it for each entry in element solution
    auto elemSol = elementSolution(element, sol, fvGeometry.gridGeometry());
    const auto centerSol = evalSolution(element, elemGeom, fvGeometry.gridGeometry(), elemSol, pos);
    for (unsigned int i = 0; i < fvGeometry.numScv(); ++i)
        elemSol[i] = centerSol;
 
    // Update volume variables with the interpolated solution. Note that this standard
    // implementation only works for element-wise constant parameters as we simply use
    // the first element scv for the vol var update. For heterogeneities within the element
    // or more complex models (e.g. 2p with interface solver) a corresponding overload
    // of this function has to be provided!
    volVars.update(elemSol, problem, element, *scvs(fvGeometry).begin());
}
} // end namespace FacetCoupling
 
template< class MDTraits,
          class CouplingMapper,
          std::size_t bulkDomainId = 0,
          std::size_t lowDimDomainId = 1,
          DiscretizationMethod bulkDM = GetPropType<typename MDTraits::template SubDomain<bulkDomainId>::TypeTag, Properties::GridGeometry>::discMethod >
class FacetCouplingManager;
 
template< class MDTraits,
          class CouplingMapper,
          std::size_t bulkDomainId = 0,
          std::size_t facetDomainId = 1,
          std::size_t edgeDomainId = 2 >
class FacetCouplingThreeDomainManager
: public FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, facetDomainId>
, public FacetCouplingManager<MDTraits, CouplingMapper, facetDomainId, edgeDomainId>
{
    using BulkFacetManager = FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, facetDomainId>;
    using FacetEdgeManager = FacetCouplingManager<MDTraits, CouplingMapper, facetDomainId, edgeDomainId>;
 
    // convenience aliases and instances of the domain ids
    using BulkIdType = typename MDTraits::template SubDomain<bulkDomainId>::Index;
    using FacetIdType = typename MDTraits::template SubDomain<facetDomainId>::Index;
    using EdgeIdType = typename MDTraits::template SubDomain<edgeDomainId>::Index;
    static constexpr auto bulkId = BulkIdType();
    static constexpr auto facetId = FacetIdType();
    static constexpr auto edgeId = EdgeIdType();
 
    // the sub-domain type tags
    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
    // further types specific to the sub-problems
    template<std::size_t id> using LocalResidual = GetPropType<SubDomainTypeTag<id>, Properties::LocalResidual>;
    template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
    template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
 
    template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
    template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
    template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::SubControlVolumeFace;
    template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
    template<std::size_t id> using GridIndexType = typename IndexTraits<GridView<id>>::GridIndex;
    template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
 
    template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
    template<std::size_t id> using ElementVolumeVariables = typename GridVariables<id>::GridVolumeVariables::LocalView;
    template<std::size_t id> using ElementFluxVariablesCache = typename GridVariables<id>::GridFluxVariablesCache::LocalView;
 
    // helper function to check if a domain uses mpfa
    template<std::size_t id>
    static constexpr bool usesMpfa(Dune::index_constant<id> domainId)
    { return GridGeometry<domainId>::discMethod == DiscretizationMethod::ccmpfa; }
 
public:
    template<std::size_t i, std::size_t j>
    using CouplingStencilType = typename std::conditional< (j == edgeDomainId),
                                                            typename FacetEdgeManager::template CouplingStencilType<i, j>,
                                                            typename BulkFacetManager::template CouplingStencilType<i, j> >::type;
 
    using SolutionVector = typename MDTraits::SolutionVector;
 
    void init(std::shared_ptr< Problem<bulkId> > bulkProblem,
              std::shared_ptr< Problem<facetId> > facetProblem,
              std::shared_ptr< Problem<edgeId> > edgeProblem,
              std::shared_ptr< CouplingMapper > couplingMapper,
              const SolutionVector& curSol)
    {
        BulkFacetManager::init(bulkProblem, facetProblem, couplingMapper, curSol);
        FacetEdgeManager::init(facetProblem, edgeProblem, couplingMapper, curSol);
    }
 
    using BulkFacetManager::couplingStencil;
    using FacetEdgeManager::couplingStencil;
 
    using BulkFacetManager::isCoupled;
    using FacetEdgeManager::isCoupled;
 
    using BulkFacetManager::isOnInteriorBoundary;
    using FacetEdgeManager::isOnInteriorBoundary;
 
    using BulkFacetManager::getLowDimVolVars;
    using FacetEdgeManager::getLowDimVolVars;
 
    using BulkFacetManager::getLowDimElement;
    using FacetEdgeManager::getLowDimElement;
 
    using BulkFacetManager::getLowDimElementIndex;
    using FacetEdgeManager::getLowDimElementIndex;
 
    using BulkFacetManager::evalSourcesFromBulk;
    using FacetEdgeManager::evalSourcesFromBulk;
 
    using BulkFacetManager::evalCouplingResidual;
    using FacetEdgeManager::evalCouplingResidual;
 
    using BulkFacetManager::bindCouplingContext;
    using FacetEdgeManager::bindCouplingContext;
 
    using BulkFacetManager::updateCouplingContext;
    using FacetEdgeManager::updateCouplingContext;
 
    using BulkFacetManager::updateCoupledVariables;
    using FacetEdgeManager::updateCoupledVariables;
 
    using BulkFacetManager::extendJacobianPattern;
    using FacetEdgeManager::extendJacobianPattern;
 
    using BulkFacetManager::evalAdditionalDomainDerivatives;
    using FacetEdgeManager::evalAdditionalDomainDerivatives;
 
    // extension of the jacobian pattern for the facet domain only occurs
    // within the bulk-facet coupling & for mpfa being used in the bulk domain.
    template<class JacobianPattern>
    void extendJacobianPattern(FacetIdType, JacobianPattern& pattern) const
    { BulkFacetManager::extendJacobianPattern(facetId, pattern); }
 
    template<class FacetLocalAssembler, class JacobianMatrixDiagBlock, class GridVariables>
    void evalAdditionalDomainDerivatives(FacetIdType,
                                         const FacetLocalAssembler& facetLocalAssembler,
                                         const typename FacetLocalAssembler::LocalResidual::ElementResidualVector& origResiduals,
                                         JacobianMatrixDiagBlock& A,
                                         GridVariables& gridVariables)
    { BulkFacetManager::evalAdditionalDomainDerivatives(facetId, facetLocalAssembler, origResiduals, A, gridVariables); }
 
    // Overload required for compatibility with mpfa
    template<class GetTensor, bool mpfa = usesMpfa(bulkId), std::enable_if_t<mpfa, int> = 0>
    auto getLowDimTensor(const Element<bulkId>& element,
                         const SubControlVolumeFace<bulkId>& scvf,
                         const GetTensor& getT) const
    { return BulkFacetManager::getLowDimTensor(element, scvf, getT); }
 
    // Overload required for compatibility with mpfa
    template<class GetTensor, bool mpfa = usesMpfa(facetId), std::enable_if_t<mpfa, int> = 0>
    auto getLowDimTensor(const Element<facetId>& element,
                         const SubControlVolumeFace<facetId>& scvf,
                         const GetTensor& getT) const
    { return FacetEdgeManager::getLowDimTensor(element, scvf, getT); }
 
    const CouplingStencilType<bulkId, edgeId>& couplingStencil(BulkIdType domainI,
                                                               const Element<bulkId>& element,
                                                               EdgeIdType domainJ) const
    { return FacetEdgeManager::getEmptyStencil(edgeId); }
 
    const CouplingStencilType<edgeId, bulkId>& couplingStencil(EdgeIdType domainI,
                                                               const Element<edgeId>& element,
                                                               BulkIdType domainJ) const
    { return BulkFacetManager::getEmptyStencil(bulkId); }
 
    void updateSolution(const SolutionVector& sol)
    {
        BulkFacetManager::updateSolution(sol);
        FacetEdgeManager::updateSolution(sol);
    }
 
    template<std::size_t i,
             std::size_t j,
             class LocalAssembler,
             std::enable_if_t<((i==bulkId && j==edgeId) || ((i==edgeId && j==bulkId))), int> = 0>
    typename LocalResidual<i>::ElementResidualVector
    evalCouplingResidual(Dune::index_constant<i> domainI,
                         const LocalAssembler& localAssembler,
                         Dune::index_constant<j> domainJ,
                         GridIndexType<j> dofIdxGlobalJ)
    {
        typename LocalResidual<i>::ElementResidualVector res(1);
        res = 0.0;
        return res;
    }
 
    template< class Assembler >
    void bindCouplingContext(FacetIdType, const Element<facetId>& element, const Assembler& assembler)
    {
        BulkFacetManager::bindCouplingContext(facetId, element, assembler);
        FacetEdgeManager::bindCouplingContext(facetId, element, assembler);
    }
 
    template< class FacetLocalAssembler >
    void updateCouplingContext(FacetIdType domainI,
                               const FacetLocalAssembler& facetLocalAssembler,
                               FacetIdType domainJ,
                               GridIndexType<facetId> dofIdxGlobalJ,
                               const PrimaryVariables<facetId>& priVarsJ,
                               unsigned int pvIdxJ)
    {
        BulkFacetManager::updateCouplingContext(domainI, facetLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
        FacetEdgeManager::updateCouplingContext(domainI, facetLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
    }
 
    template<std::size_t i,
             std::size_t j,
             class LocalAssembler,
             std::enable_if_t<((i==bulkId && j==edgeId) || (i==edgeId && j==bulkId)), int> = 0>
    void updateCouplingContext(Dune::index_constant<i> domainI,
                               const LocalAssembler& localAssembler,
                               Dune::index_constant<j> domainJ,
                               GridIndexType<j> dofIdxGlobalJ,
                               const PrimaryVariables<j>& priVarsJ,
                               unsigned int pvIdxJ)
    { /*do nothing here*/ }
 
    template< class FacetLocalAssembler, class UpdatableElementVolVars, class UpdatableFluxVarCache>
    void updateCoupledVariables(FacetIdType domainI,
                                const FacetLocalAssembler& facetLocalAssembler,
                                UpdatableElementVolVars& elemVolVars,
                                UpdatableFluxVarCache& elemFluxVarsCache)
    {
        BulkFacetManager::updateCoupledVariables(domainI, facetLocalAssembler, elemVolVars, elemFluxVarsCache);
        FacetEdgeManager::updateCoupledVariables(domainI, facetLocalAssembler, elemVolVars, elemFluxVarsCache);
    }
 
    template<std::size_t id, std::enable_if_t<(id == bulkId || id == facetId), int> = 0>
    const Problem<id>& problem() const { return BulkFacetManager::template problem<id>(); }
 
    template<std::size_t id, std::enable_if_t<(id == bulkId || id == facetId), int> = 0>
    Problem<id>& problem() { return BulkFacetManager::template problem<id>(); }
 
    template<std::size_t id, std::enable_if_t<(id == edgeId), int> = 0>
    const Problem<id>& problem() const { return FacetEdgeManager::template problem<id>(); }
 
    template<std::size_t id, std::enable_if_t<(id == edgeId), int> = 0>
    Problem<id>& problem() { return FacetEdgeManager::template problem<id>(); }
};
 
} // end namespace Dumux
 
// Here, we have to include all available implementations
#include <dumux/multidomain/facet/box/couplingmanager.hh>
#include <dumux/multidomain/facet/cellcentered/tpfa/couplingmanager.hh>
#include <dumux/multidomain/facet/cellcentered/mpfa/couplingmanager.hh>
 
#endif