diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh
index 938ae4915987f8fc2b1da60d67b9466c2125c6ae..cb3b2b57ade83aa1c4952885cf3f9b9d1ad79bec 100644
--- a/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh
+++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh
@@ -82,13 +82,6 @@ constexpr auto coupledDomains(Dune::index_constant<i> domainI)
return std::make_tuple(freeFlowMomentumIndex, freeFlowMassIndex);
}
-template<std::size_t i, std::size_t j>
-struct LocalIndices
-{
- static constexpr auto domainI = Dune::index_constant<i>();
- static constexpr auto domainJ = Dune::index_constant<j>();
-};
-
template<std::size_t i, std::size_t j>
constexpr auto globalToLocalDomainIndices(Dune::index_constant<i>, Dune::index_constant<j>)
{
@@ -96,20 +89,20 @@ constexpr auto globalToLocalDomainIndices(Dune::index_constant<i>, Dune::index_c
static_assert(i != j);
if constexpr (i < j)
- return LocalIndices<0, 1>();
+ return std::pair<Dune::index_constant<0>, Dune::index_constant<1>>{};
else
- return LocalIndices<1, 0>();
+ return std::pair<Dune::index_constant<1>, Dune::index_constant<0>>{};
}
-struct CouplingManagerMaps
+struct CouplingMaps
{
- static constexpr auto makeCouplingManagerMap()
+ static constexpr auto managerMap()
{
return FreeFlowPorousMediumDetail::makeCouplingManagerMap();
}
template<std::size_t i, std::size_t j>
- static constexpr auto globalToLocalDomainIndices(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
+ static constexpr auto globalToLocal(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
{
return FreeFlowPorousMediumDetail::globalToLocalDomainIndices(domainI, domainJ);
}
@@ -158,7 +151,7 @@ template<class MDTraits>
class FreeFlowPorousMediumCouplingManager
: public MultiBinaryCouplingManager<
MDTraits,
- FreeFlowPorousMediumDetail::CouplingManagerMaps,
+ FreeFlowPorousMediumDetail::CouplingMaps,
typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowCouplingManager,
typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMomentumPorousMediumCouplingManager,
typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMassPorousMediumCouplingManager
@@ -166,7 +159,7 @@ class FreeFlowPorousMediumCouplingManager
{
using ParentType = MultiBinaryCouplingManager<
MDTraits,
- FreeFlowPorousMediumDetail::CouplingManagerMaps,
+ FreeFlowPorousMediumDetail::CouplingMaps,
typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowCouplingManager,
typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMomentumPorousMediumCouplingManager,
typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMassPorousMediumCouplingManager
@@ -202,6 +195,8 @@ public:
static constexpr auto porousMediumIndex = FreeFlowPorousMediumDetail::porousMediumIndex;
public:
+ using ParentType::ParentType;
+
template<class GridVarsTuple>
void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem,
std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
@@ -247,9 +242,9 @@ public:
{
static_assert(domainI != freeFlowMomentumIndex && domainJ != freeFlowMomentumIndex);
- auto& couplingContext = this->subCouplingManager(domainI, domainJ).couplingContext(
- ParentType::globalToLocalDomainIndices(domainI, domainJ).domainI, fvGeometry, scvf
- );
+ const auto& couplingContext = this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
+ return cm.couplingContext(ii, fvGeometry, scvf);
+ });
const auto& freeFlowElement = [&]
{
@@ -358,8 +353,9 @@ public:
const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
const bool considerPreviousTimeStep = false) const
{
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex)
- .density(element, fvGeometry, scvf, considerPreviousTimeStep);
+ return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).density(
+ element, fvGeometry, scvf, considerPreviousTimeStep
+ );
}
auto insideAndOutsideDensity(const Element<freeFlowMomentumIndex>& element,
@@ -407,43 +403,12 @@ public:
);
}
-
- template<std::size_t i>
- const auto& problem(Dune::index_constant<i> domainI) const
- {
- if constexpr (domainI == freeFlowMomentumIndex)
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).problem(
- ParentType::globalToLocalDomainIndices(freeFlowMomentumIndex, freeFlowMassIndex).domainI
- );
- else if constexpr (domainI == freeFlowMassIndex)
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).problem(
- ParentType::globalToLocalDomainIndices(freeFlowMomentumIndex, freeFlowMassIndex).domainJ
- );
- else if constexpr (domainI == porousMediumIndex)
- return this->subCouplingManager(porousMediumIndex, freeFlowMassIndex).problem(
- ParentType::globalToLocalDomainIndices(porousMediumIndex, freeFlowMassIndex).domainI
- );
- else
- DUNE_THROW(Dune::InvalidStateException, "Wrong index");
- }
-
template<std::size_t i>
- auto& problem(Dune::index_constant<i> domainI)
+ const Problem<i>& problem(Dune::index_constant<i> domainI) const
{
- if constexpr (domainI == freeFlowMomentumIndex)
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).problem(
- ParentType::globalToLocalDomainIndices(freeFlowMomentumIndex, freeFlowMassIndex).domainI
- );
- else if constexpr (domainI == freeFlowMassIndex)
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).problem(
- ParentType::globalToLocalDomainIndices(freeFlowMomentumIndex, freeFlowMassIndex).domainJ
- );
- else if constexpr (domainI == porousMediumIndex)
- return this->subCouplingManager(porousMediumIndex, freeFlowMassIndex).problem(
- ParentType::globalToLocalDomainIndices(porousMediumIndex, freeFlowMassIndex).domainI
- );
- else
- DUNE_THROW(Dune::InvalidStateException, "Wrong index");
+ return this->subApply(domainI, [&](const auto& cm, auto&& ii) -> const auto& {
+ return cm.problem(ii);
+ });
}
template<std::size_t i, std::size_t j>
@@ -451,9 +416,9 @@ public:
Dune::index_constant<j> domainJ,
const SubControlVolumeFace<i>& scvf) const
{
- return this->subCouplingManager(domainI, domainJ).isCoupled(
- ParentType::globalToLocalDomainIndices(domainI, domainJ).domainI, scvf
- );
+ return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
+ return cm.isCoupled(ii, scvf);
+ });
}
/*!
@@ -466,9 +431,9 @@ public:
Dune::index_constant<j> domainJ,
const SubControlVolume<i>& scv) const
{
- return this->subCouplingManager(domainI, domainJ).isCoupled(
- ParentType::globalToLocalDomainIndices(domainI, domainJ).domainI, scv
- );
+ return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
+ return cm.isCoupled(ii, scv);
+ });
}
/*!
@@ -478,9 +443,9 @@ public:
Dune::index_constant<porousMediumIndex> domainJ,
const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
{
- return this->subCouplingManager(domainI, domainJ).isCoupledLateralScvf(
- ParentType::globalToLocalDomainIndices(domainI, domainJ).domainI, scvf
- );
+ return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
+ return cm.isCoupledLateralScvf(ii, scvf);
+ });
}
@@ -501,11 +466,9 @@ public:
Dune::index_constant<j> domainJ) const
{
static_assert(freeFlowMomentumIndex != j);
- return this->subCouplingManager(domainI, domainJ).couplingStencil(
- ParentType::globalToLocalDomainIndices(domainI, domainJ).domainI,
- elementI, scvI,
- ParentType::globalToLocalDomainIndices(domainI, domainJ).domainJ
- );
+ return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
+ return cm.couplingStencil(ii, elementI, scvI, jj);
+ });
}
};
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh
index 09b317561d8e958c01e1701f4bd0aa0a5db2bb01..e8156b97747aff358f6f33891fcb029db0f2380d 100644
--- a/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh
@@ -117,11 +117,9 @@ public:
std::shared_ptr<Problem<porousMediumIndex>> darcyProblem,
SolutionVectorStorage& curSol)
{
+ this->setSubProblems(std::make_tuple(freeFlowMassProblem, darcyProblem));
this->attachSolution(curSol);
- freeFlowMassProblemPtr_ = &(*freeFlowMassProblem);
- porousMediumProblemPtr_ = &(*darcyProblem);
-
couplingMapper_.update(*this);
}
@@ -164,9 +162,9 @@ public:
{
if (c.dofIdx == dofIdxGlobalJ)
{
- const auto elemSol = elementSolution(c.element, this->curSol(domainJ), problem(domainJ).gridGeometry());
+ const auto elemSol = elementSolution(c.element, this->curSol(domainJ), this->problem(domainJ).gridGeometry());
const auto& scv = *scvs(c.fvGeometry).begin();
- c.volVars.update(elemSol, problem(domainJ), c.element, scv);
+ c.volVars.update(elemSol, this->problem(domainJ), c.element, scv);
}
}
}
@@ -209,7 +207,7 @@ public:
const Element<i>& element,
Dune::index_constant<j> domainJ) const
{
- const auto eIdx = problem(domainI).gridGeometry().elementMapper().index(element);
+ const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
return couplingMapper_.couplingStencil(domainI, eIdx, domainJ);
}
@@ -222,24 +220,6 @@ public:
bool isCoupled(Dune::index_constant<i> domainI, const SubControlVolumeFace<i>& scvf) const
{ return couplingMapper_.isCoupled(domainI, scvf); }
- template<std::size_t i>
- Problem<i>& problem(Dune::index_constant<i> domainI)
- {
- if constexpr (i == freeFlowMassIndex)
- return *freeFlowMassProblemPtr_;
- else
- return *porousMediumProblemPtr_;
- }
-
- template<std::size_t i>
- const Problem<i>& problem(Dune::index_constant<i> domainI) const
- {
- if constexpr (i == freeFlowMassIndex)
- return *freeFlowMassProblemPtr_;
- else
- return *porousMediumProblemPtr_;
- }
-
private:
/*!
* \brief Returns whether a given scvf is coupled to the other domain
@@ -268,7 +248,7 @@ private:
template<std::size_t i>
void bindCouplingContext_(Dune::index_constant<i> domainI, const Element<i>& element) const
{
- const auto fvGeometry = localView(problem(domainI).gridGeometry()).bindElement(element);
+ const auto fvGeometry = localView(this->problem(domainI).gridGeometry()).bindElement(element);
bindCouplingContext_(domainI, fvGeometry);
}
@@ -281,7 +261,7 @@ private:
auto& context = std::get<domainI>(couplingContext_);
context.clear();
- const auto eIdx = problem(domainI).gridGeometry().elementMapper().index(fvGeometry.element());
+ const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(fvGeometry.element());
// do nothing if the element is not coupled to the other domain
if (!isCoupledElement_(domainI, eIdx))
@@ -295,7 +275,7 @@ private:
{
const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf);
constexpr auto domainJ = Dune::index_constant<1-domainI>();
- const auto& otherGridGeometry = problem(domainJ).gridGeometry();
+ const auto& otherGridGeometry = this->problem(domainJ).gridGeometry();
const auto& otherElement = otherGridGeometry.element(otherElementIdx);
auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement);
@@ -316,9 +296,6 @@ private:
mutable std::tuple<std::vector<FreeFlowMassCouplingContext>, std::vector<PorousMediumCouplingContext>> couplingContext_;
mutable std::array<std::size_t, 2> couplingContextBoundForElement_;
- Problem<freeFlowMassIndex>* freeFlowMassProblemPtr_ = nullptr;
- Problem<porousMediumIndex>* porousMediumProblemPtr_ = nullptr;
-
CouplingMapper couplingMapper_;
};
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh
index c8193ac26430470d82c1e7ff294bb1575290b0c2..f957e9f65f15d9c92dac5e867681b8f63f7d14f3 100644
--- a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh
@@ -116,11 +116,9 @@ public:
std::shared_ptr<Problem<porousMediumIndex>> porousMediumProblem,
SolutionVectorStorage& curSol)
{
+ this->setSubProblems(std::make_tuple(freeFlowMomentumProblem, porousMediumProblem));
this->attachSolution(curSol);
- freeFlowMomentumProblemPtr_ = &(*freeFlowMomentumProblem);
- porousMediumProblemPtr_ = &(*porousMediumProblem);
-
couplingMapper_.update(*this);
}
@@ -187,7 +185,7 @@ public:
const auto& ggJ = c.fvGeometry.gridGeometry();
const auto& scv = *scvs(c.fvGeometry).begin();
const auto elemSol = elementSolution(c.element, this->curSol(domainJ), ggJ);
- c.volVars.update(elemSol, problem(domainJ), c.element, scv);
+ c.volVars.update(elemSol, this->problem(domainJ), c.element, scv);
}
}
}
@@ -230,7 +228,7 @@ public:
const Element<porousMediumIndex>& element,
Dune::index_constant<freeFlowMomentumIndex> domainJ) const
{
- const auto eIdx = problem(domainI).gridGeometry().elementMapper().index(element);
+ const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
return couplingMapper_.couplingStencil(domainI, eIdx, domainJ);
}
@@ -315,24 +313,6 @@ public:
return velocity;
}
- template<std::size_t i>
- Problem<i>& problem(Dune::index_constant<i> domainI)
- {
- if constexpr (i == freeFlowMomentumIndex)
- return *freeFlowMomentumProblemPtr_;
- else
- return *porousMediumProblemPtr_;
- }
-
- template<std::size_t i>
- const Problem<i>& problem(Dune::index_constant<i> domainI) const
- {
- if constexpr (i == freeFlowMomentumIndex)
- return *freeFlowMomentumProblemPtr_;
- else
- return *porousMediumProblemPtr_;
- }
-
private:
//! Return the volume variables of domain i for a given element and scv
template<std::size_t i>
@@ -342,9 +322,9 @@ private:
{
VolumeVariables<i> volVars;
const auto elemSol = elementSolution(
- element, this->curSol(domainI), problem(domainI).gridGeometry()
+ element, this->curSol(domainI), this->problem(domainI).gridGeometry()
);
- volVars.update(elemSol, problem(domainI), element, scv);
+ volVars.update(elemSol, this->problem(domainI), element, scv);
return volVars;
}
@@ -361,7 +341,7 @@ private:
template<std::size_t i>
void bindCouplingContext_(Dune::index_constant<i> domainI, const Element<i>& element) const
{
- const auto fvGeometry = localView(problem(domainI).gridGeometry()).bindElement(element);
+ const auto fvGeometry = localView(this->problem(domainI).gridGeometry()).bindElement(element);
bindCouplingContext_(domainI, fvGeometry);
}
@@ -374,7 +354,7 @@ private:
auto& context = std::get<domainI>(couplingContext_);
context.clear();
- const auto eIdx = problem(domainI).gridGeometry().elementMapper().index(fvGeometry.element());
+ const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(fvGeometry.element());
// do nothing if the element is not coupled to the other domain
if (!isCoupledElement_(domainI, eIdx))
@@ -390,7 +370,7 @@ private:
{
const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf);
constexpr auto domainJ = Dune::index_constant<1-i>();
- const auto& otherGridGeometry = problem(domainJ).gridGeometry();
+ const auto& otherGridGeometry = this->problem(domainJ).gridGeometry();
const auto& otherElement = otherGridGeometry.element(otherElementIdx);
auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement);
@@ -402,7 +382,7 @@ private:
scvf.index(),
couplingMapper_.flipScvfIndex(domainI, scvf),
otherElementIdx,
- problem(domainJ).spatialParams().gravity(scvf.center())
+ this->problem(domainJ).spatialParams().gravity(scvf.center())
});
}
@@ -410,7 +390,7 @@ private:
{
const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf);
constexpr auto domainJ = Dune::index_constant<1-i>();
- const auto& otherGridGeometry = problem(domainJ).gridGeometry();
+ const auto& otherGridGeometry = this->problem(domainJ).gridGeometry();
const auto& otherElement = otherGridGeometry.element(otherElementIdx);
auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement);
@@ -434,9 +414,6 @@ private:
mutable std::tuple<std::vector<FreeFlowMomentumCouplingContext>, std::vector<PorousMediumCouplingContext>> couplingContext_;
mutable std::array<std::size_t, 2> couplingContextBoundForElement_;
- Problem<freeFlowMomentumIndex>* freeFlowMomentumProblemPtr_ = nullptr;
- Problem<porousMediumIndex>* porousMediumProblemPtr_ = nullptr;
-
CouplingMapper couplingMapper_;
};
diff --git a/dumux/multidomain/multibinarycouplingmanager.hh b/dumux/multidomain/multibinarycouplingmanager.hh
index 082c69018766909e815e2e7204b8880fe28bf655..cd335ffb5bc48361b169ae87d43f27e9512f2df7 100644
--- a/dumux/multidomain/multibinarycouplingmanager.hh
+++ b/dumux/multidomain/multibinarycouplingmanager.hh
@@ -43,17 +43,21 @@ struct HasIndex<i, std::tuple<Indices...>>
: std::disjunction<std::is_same<Dune::index_constant<i>, Indices>...>
{};
-template<class Map, std::size_t i, std::size_t j>
-constexpr bool isCoupled(Dune::index_constant<i> I, Dune::index_constant<j>)
-{ return HasIndex<j, std::decay_t<decltype(Map::coupledDomains(I))>>::value; }
-
} // end namespace Detail
/*!
* \ingroup MultiDomain
* \brief Coupling manager that combines an arbitrary number of binary coupling manager (coupling two domains each)
+ * \tparam MDTraits the multidomain traits
+ * \tparam CouplingMap a coupling policy class
+ * \tparam CouplingMgrs the binary sub-coupling manager types
+ *
+ * The coupling policy has to provide the interfaces
+ * - CouplingMap::coupledDomains(i): returns a tuple of Dune::index_constants with the coupled domains
+ * - CouplingMap::globalToLocal(i, j): maps the indices i, j to the local index pair of the responsible sub coupling manager
+ * - CouplingMap::managerMap(): returns a two-dimensional array mapping two indices to the coupling manager index
*/
-template<class MDTraits, class CouplingManagerMaps, class ...CouplingMgrs>
+template<class MDTraits, class CouplingMap, class ...CouplingMgrs>
class MultiBinaryCouplingManager
{
// the sub domain type tags
@@ -78,12 +82,35 @@ class MultiBinaryCouplingManager
using SubSolutionVector = std::decay_t<decltype(std::declval<typename MDTraits::SolutionVector>()[Dune::index_constant<id>()])>;
using SolutionVectors = typename MDTraits::template TupleOfSharedPtr<SubSolutionVector>;
-public:
+ static constexpr auto couplingManagerMap_ = CouplingMap::managerMap();
+
+ //! Returns the local domain indices used within the binary sub coupling managers.
+ template<std::size_t i, std::size_t j>
+ static constexpr auto globalToLocal_(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
+ {
+ static_assert(i <= MDTraits::numSubDomains && j <= MDTraits::numSubDomains);
+ return CouplingMap::globalToLocal(domainI, domainJ);
+ }
- static constexpr auto couplingManagerMap = CouplingManagerMaps::makeCouplingManagerMap();
+ //! If two domain are coupled
+ template<class Map, std::size_t i, std::size_t j>
+ static constexpr bool isCoupled_(Dune::index_constant<i> domainI, Dune::index_constant<j>)
+ { return Detail::HasIndex<j, std::decay_t<decltype(Map::coupledDomains(domainI))>>::value; }
+ //! Returns the coupling manager index for a given domain combination
template<std::size_t i, std::size_t j>
- using SubCouplingManager = SubCouplingManagerT<couplingManagerMap[i][j]>;
+ static constexpr auto subCouplingManagerIndex_(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
+ {
+ static_assert(
+ isCoupled_<CouplingMap>(domainI, domainJ),
+ "Sub-coupling manager only exists for coupled domains."
+ );
+ return couplingManagerMap_[i][j];
+ }
+
+public:
+ template<std::size_t i, std::size_t j>
+ using SubCouplingManager = SubCouplingManagerT<couplingManagerMap_[i][j]>;
MultiBinaryCouplingManager()
{
@@ -99,31 +126,12 @@ public:
});
}
- //! Returns the local domain indices used within the binary sub coupling managers.
- template<std::size_t i, std::size_t j>
- static constexpr auto globalToLocalDomainIndices(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
- {
- static_assert(i <= MDTraits::numSubDomains && j <= MDTraits::numSubDomains);
- return CouplingManagerMaps::globalToLocalDomainIndices(domainI, domainJ);
- }
-
- //! Returns the coupling manager index for a given domain combination
- template<std::size_t i, std::size_t j>
- static constexpr auto getCouplingManagerIndex(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
- {
- static_assert(
- Detail::template isCoupled<CouplingManagerMaps>(domainI, domainJ),
- "Sub-coupling manager only exists for coupled domains."
- );
- return couplingManagerMap[i][j];
- }
-
//! return the binary sub-coupling manager
template<std::size_t i, std::size_t j>
auto& subCouplingManager(Dune::index_constant<i> domainI,
Dune::index_constant<j> domainJ)
{
- constexpr auto idx = getCouplingManagerIndex(domainI, domainJ);
+ constexpr auto idx = subCouplingManagerIndex_(domainI, domainJ);
return *std::get<idx>(couplingManagers_);
}
@@ -132,10 +140,52 @@ public:
const auto& subCouplingManager(Dune::index_constant<i> domainI,
Dune::index_constant<j> domainJ) const
{
- constexpr auto idx = getCouplingManagerIndex(domainI, domainJ);
+ constexpr auto idx = subCouplingManagerIndex_(domainI, domainJ);
return *std::get<idx>(couplingManagers_);
}
+ //! apply a function to the domainI-domainJ sub coupling manager using its local indices
+ template<std::size_t i, std::size_t j, class Apply>
+ decltype(auto) subApply(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ Apply&& apply)
+ {
+ constexpr auto localIndices = globalToLocal_(domainI, domainJ);
+ return apply(subCouplingManager(domainI, domainJ), localIndices.first, localIndices.second);
+ }
+
+ //! apply a function to the domainI-domainJ sub coupling manager using its local indices
+ template<std::size_t i, std::size_t j, class Apply>
+ decltype(auto) subApply(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const Apply& apply) const
+ {
+ constexpr auto localIndices = globalToLocal_(domainI, domainJ);
+ return apply(subCouplingManager(domainI, domainJ), localIndices.first, localIndices.second);
+ }
+
+ //! apply a function to a sub coupling manager containing this domain
+ template<std::size_t i, class Apply>
+ decltype(auto) subApply(Dune::index_constant<i> domainI, Apply&& apply)
+ {
+ constexpr auto dm = CouplingMap::coupledDomains(domainI);
+ static_assert(std::tuple_size_v<std::decay_t<decltype(dm)>> != 0, "subApply on uncoupled domain!");
+ constexpr auto domainJ = std::get<0>(dm);
+ constexpr auto localIndices = globalToLocal_(domainI, domainJ);
+ return apply(subCouplingManager(domainI, domainJ), localIndices.first);
+ }
+
+ //! apply a function to a sub coupling manager containing this domain
+ template<std::size_t i, class Apply>
+ decltype(auto) subApply(Dune::index_constant<i> domainI, const Apply& apply) const
+ {
+ constexpr auto dm = CouplingMap::coupledDomains(domainI);
+ static_assert(std::tuple_size_v<std::decay_t<decltype(dm)>> != 0, "subApply on uncoupled domain!");
+ constexpr auto domainJ = std::get<0>(dm);
+ constexpr auto localIndices = globalToLocal_(domainI, domainJ);
+ return apply(subCouplingManager(domainI, domainJ), localIndices.first);
+ }
+
//! Update the solution vector before assembly
void updateSolution(const typename MDTraits::SolutionVector& curSol)
{
@@ -167,11 +217,10 @@ public:
Dune::index_constant<j> domainJ) const
{
// if the domains are coupled according to the map, forward to sub-coupling manager
- if constexpr (Detail::template isCoupled<CouplingManagerMaps>(domainI, domainJ))
- return subCouplingManager(domainI, domainJ).couplingStencil(
- globalToLocalDomainIndices(domainI, domainJ).domainI, entity,
- globalToLocalDomainIndices(domainI, domainJ).domainJ
- );
+ if constexpr (isCoupled_<CouplingMap>(domainI, domainJ))
+ return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
+ return cm.couplingStencil(ii, entity, jj);
+ });
else
return emptyStencil_;
}
@@ -185,13 +234,10 @@ public:
Dune::index_constant<j> domainJ,
std::size_t dofIdxGlobalJ) const
{
- if constexpr (i == j || Detail::template isCoupled<CouplingManagerMaps>(domainI, domainJ))
- return subCouplingManager(domainI, domainJ).evalCouplingResidual(
- globalToLocalDomainIndices(domainI, domainJ).domainI,
- scvfI, localAssemblerI,
- globalToLocalDomainIndices(domainI, domainJ).domainJ,
- dofIdxGlobalJ
- );
+ if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
+ return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
+ return cm.evalCouplingResidual(ii, scvfI, localAssemblerI, jj, dofIdxGlobalJ);
+ });
else
{
DUNE_THROW(Dune::InvalidStateException,
@@ -209,13 +255,10 @@ public:
Dune::index_constant<j> domainJ,
std::size_t dofIdxGlobalJ) const
{
- if constexpr (i == j || Detail::template isCoupled<CouplingManagerMaps>(domainI, domainJ))
- return subCouplingManager(domainI, domainJ).evalCouplingResidual(
- globalToLocalDomainIndices(domainI, domainJ).domainI,
- localAssemblerI,
- globalToLocalDomainIndices(domainI, domainJ).domainJ,
- dofIdxGlobalJ
- );
+ if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
+ return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
+ return cm.evalCouplingResidual(ii, localAssemblerI, jj, dofIdxGlobalJ);
+ });
else
{
DUNE_THROW(Dune::InvalidStateException,
@@ -234,14 +277,10 @@ public:
Dune::index_constant<j> domainJ,
std::size_t dofIdxGlobalJ) const
{
- if constexpr (i == j || Detail::template isCoupled<CouplingManagerMaps>(domainI, domainJ))
- return subCouplingManager(domainI, domainJ).evalCouplingResidual(
- globalToLocalDomainIndices(domainI, domainJ).domainI,
- localAssemblerI,
- scvI,
- globalToLocalDomainIndices(domainI, domainJ).domainJ,
- dofIdxGlobalJ
- );
+ if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
+ return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
+ return cm.evalCouplingResidual(ii, localAssemblerI, scvI, jj, dofIdxGlobalJ);
+ });
else
{
DUNE_THROW(Dune::InvalidStateException,
@@ -265,32 +304,21 @@ public:
// only one other manager needs to take action if i != j
if constexpr (i != j)
{
- if constexpr (Detail::template isCoupled<CouplingManagerMaps>(domainI, domainJ))
- subCouplingManager(domainI, domainJ).updateCouplingContext(
- globalToLocalDomainIndices(domainI, domainJ).domainI,
- localAssemblerI,
- globalToLocalDomainIndices(domainI, domainJ).domainJ,
- dofIdxGlobalJ,
- priVars,
- pvIdxJ
- );
+ if constexpr (isCoupled_<CouplingMap>(domainI, domainJ))
+ subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
+ cm.updateCouplingContext(ii, localAssemblerI, jj, dofIdxGlobalJ, priVars, pvIdxJ);
+ });
}
else
{
// for i == j, we need to call all relevant managers where domainI is involved and make sure that the
- // context is updated with respect to its own domain
+ // context is updated with respect to its own domain (I-I coupling context)
using namespace Dune::Hybrid;
- forEach(CouplingManagerMaps::coupledDomains(domainI), [&](const auto domainJ)
- {
+ forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
static_assert(domainI != domainJ);
- subCouplingManager(domainI, domainJ).updateCouplingContext(
- globalToLocalDomainIndices(domainI, domainJ).domainI,
- localAssemblerI,
- globalToLocalDomainIndices(domainI, domainJ).domainI,
- dofIdxGlobalJ,
- priVars,
- pvIdxJ
- );
+ subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
+ cm.updateCouplingContext(ii, localAssemblerI, ii, dofIdxGlobalJ, priVars, pvIdxJ);
+ });
});
}
}
@@ -301,11 +329,10 @@ public:
{
// for the coupling blocks
using namespace Dune::Hybrid;
- forEach(CouplingManagerMaps::coupledDomains(domainI), [&](const auto domainJ)
- {
- this->subCouplingManager(domainI, domainJ).bindCouplingContext(
- globalToLocalDomainIndices(domainI, domainJ).domainI, element, assembler
- );
+ forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
+ subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj) -> void {
+ cm.bindCouplingContext(ii, element, assembler);
+ });
});
}
@@ -317,8 +344,9 @@ public:
decltype(auto) numericEpsilon(Dune::index_constant<i> domainI,
const std::string& paramGroup) const
{
- static constexpr auto someOtherDomain = std::get<0>(CouplingManagerMaps::coupledDomains(domainI));
- return subCouplingManager(domainI, someOtherDomain).numericEpsilon(globalToLocalDomainIndices(domainI, someOtherDomain).domainI, paramGroup);
+ return subApply(domainI, [&](const auto& cm, auto&& ii) -> decltype(auto) {
+ return cm.numericEpsilon(ii, paramGroup);
+ });
}
/*!
@@ -343,14 +371,10 @@ public:
{
// for the coupling blocks
using namespace Dune::Hybrid;
- forEach(CouplingManagerMaps::coupledDomains(domainI), [&](const auto domainJ)
- {
- subCouplingManager(domainI, domainJ).updateCoupledVariables(
- globalToLocalDomainIndices(domainI, domainJ).domainI,
- localAssemblerI,
- elemVolVars,
- elemFluxVarsCache
- );
+ forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
+ subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
+ cm.updateCoupledVariables(ii, localAssemblerI, elemVolVars, elemFluxVarsCache);
+ });
});
}