Commit 004c3cd5 authored by Martin Schneider's avatar Martin Schneider
Browse files

[md][ffpm][box] Convenience functions to iterate over coupling segments

parent edfefcf4
......@@ -285,47 +285,42 @@ public:
boundDarcyElemIdx_ = darcyElementIdx;
// do nothing if the element is not coupled to the other domain
if(!couplingMapper_.pmCouplingFacetIdxMap().count(darcyElementIdx))
if(!couplingMapper_.couplingFacetIdxMap(domainI).count(darcyElementIdx))
return;
// prepare the coupling context
const auto& couplingFacets = couplingMapper_.pmCouplingFacetIdxMap().at(darcyElementIdx);
auto stokesFvGeometry = localView(this->problem(freeFlowIdx).gridGeometry());
// Iterate over facets belonging to scvfs
for(const auto& couplingFacetsScvf : couplingFacets)
for(const auto& couplingFacet : couplingFacets(domainI, couplingMapper_, darcyElementIdx))
{
for(const auto& facetIdx : couplingFacetsScvf)
{
const auto& couplingFacet = couplingMapper_.couplingFacet(facetIdx);
const auto& stokesElement = this->problem(freeFlowIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingFacet.ffEIdx);
stokesFvGeometry.bindElement(stokesElement);
const auto& stokesElement = this->problem(freeFlowIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingFacet.ffEIdx);
stokesFvGeometry.bindElement(stokesElement);
VelocityVector faceVelocity(0.0);
VelocityVector faceVelocity(0.0);
for(const auto& scvf : scvfs(stokesFvGeometry))
{
if(scvf.index() == couplingFacet.ffScvfIdx)
faceVelocity[scvf.directionIndex()] = this->curSol()[freeFlowFaceIdx][scvf.dofIndex()];
}
for(const auto& scvf : scvfs(stokesFvGeometry))
{
if(scvf.index() == couplingFacet.ffScvfIdx)
faceVelocity[scvf.directionIndex()] = this->curSol()[freeFlowFaceIdx][scvf.dofIndex()];
}
using PriVarsType = typename VolumeVariables<freeFlowCellCenterIdx>::PrimaryVariables;
const auto& cellCenterPriVars = this->curSol()[freeFlowCellCenterIdx][couplingFacet.ffEIdx];
const auto elemSol = makeElementSolutionFromCellCenterPrivars<PriVarsType>(cellCenterPriVars);
using PriVarsType = typename VolumeVariables<freeFlowCellCenterIdx>::PrimaryVariables;
const auto& cellCenterPriVars = this->curSol()[freeFlowCellCenterIdx][couplingFacet.ffEIdx];
const auto elemSol = makeElementSolutionFromCellCenterPrivars<PriVarsType>(cellCenterPriVars);
VolumeVariables<freeFlowIdx> stokesVolVars;
for(const auto& scv : scvs(stokesFvGeometry))
stokesVolVars.update(elemSol, this->problem(freeFlowIdx), stokesElement, scv);
VolumeVariables<freeFlowIdx> stokesVolVars;
for(const auto& scv : scvs(stokesFvGeometry))
stokesVolVars.update(elemSol, this->problem(freeFlowIdx), stokesElement, scv);
// add the context
darcyCouplingContext_.push_back({stokesElement, stokesFvGeometry,
couplingFacet.ffScvfIdx, couplingFacet.pmScvfIdx,
faceVelocity, stokesVolVars,
std::make_unique< std::vector<StokesCouplingContext> >(),
couplingFacet.geometry});
// add the context
darcyCouplingContext_.push_back({stokesElement, stokesFvGeometry,
couplingFacet.ffScvfIdx, couplingFacet.pmScvfIdx,
faceVelocity, stokesVolVars,
std::make_unique< std::vector<StokesCouplingContext> >(),
couplingFacet.geometry});
fillCouplingContext(stokesElement, assembler, (*darcyCouplingContext_.back().stokesContext));
}
fillCouplingContext(stokesElement, assembler, (*darcyCouplingContext_.back().stokesContext));
}
}
......@@ -831,35 +826,31 @@ private:
{
const auto stokesElementIdx = this->problem(freeFlowIdx).gridGeometry().elementMapper().index(element);
auto domainI = Dune::index_constant<freeFlowIdx>();
// do nothing if the element is not coupled to the other domain
if(!couplingMapper_.ffCouplingFacetIdxMap().count(stokesElementIdx))
if(!couplingMapper_.couplingFacetIdxMap(domainI).count(stokesElementIdx))
return;
// prepare the coupling context
const auto& couplingFacets = couplingMapper_.ffCouplingFacetIdxMap().at(stokesElementIdx);
auto darcyFvGeometry = localView(this->problem(porousMediumIdx).gridGeometry());
// Iterate over facets belonging to scvfs
for(const auto& couplingFacetsScvf : couplingFacets)
for(const auto& couplingFacet : couplingFacets(domainI, couplingMapper_, stokesElementIdx))
{
for(const auto& facetIdx : couplingFacetsScvf)
{
const auto& couplingFacet = couplingMapper_.couplingFacet(facetIdx);
const auto& darcyElement = this->problem(porousMediumIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingFacet.pmEIdx);
darcyFvGeometry.bind(darcyElement);
const auto& darcyElement = this->problem(porousMediumIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingFacet.pmEIdx);
darcyFvGeometry.bind(darcyElement);
auto darcyElemVolVars = localView(assembler.gridVariables(porousMediumIdx).curGridVolVars());
auto darcyElemFluxVarsCache = localView(assembler.gridVariables(porousMediumIdx).gridFluxVarsCache());
auto darcyElemVolVars = localView(assembler.gridVariables(porousMediumIdx).curGridVolVars());
auto darcyElemFluxVarsCache = localView(assembler.gridVariables(porousMediumIdx).gridFluxVarsCache());
darcyElemVolVars.bind(darcyElement, darcyFvGeometry, this->curSol()[porousMediumIdx]);
darcyElemFluxVarsCache.bind(darcyElement, darcyFvGeometry, darcyElemVolVars);
darcyElemVolVars.bind(darcyElement, darcyFvGeometry, this->curSol()[porousMediumIdx]);
darcyElemFluxVarsCache.bind(darcyElement, darcyFvGeometry, darcyElemVolVars);
// add the context
couplingContext.push_back({darcyElement, darcyFvGeometry, couplingFacet.pmScvfIdx, couplingFacet.ffScvfIdx,
std::make_unique< ElementVolumeVariables<porousMediumIdx> >( std::move(darcyElemVolVars)),
std::make_unique< ElementFluxVariablesCache<porousMediumIdx> >( std::move(darcyElemFluxVarsCache)),
couplingFacet.geometry});
}
// add the context
couplingContext.push_back({darcyElement, darcyFvGeometry, couplingFacet.pmScvfIdx, couplingFacet.ffScvfIdx,
std::make_unique< ElementVolumeVariables<porousMediumIdx> >( std::move(darcyElemVolVars)),
std::make_unique< ElementFluxVariablesCache<porousMediumIdx> >( std::move(darcyElemFluxVarsCache)),
couplingFacet.geometry});
}
}
......
......@@ -42,19 +42,19 @@
namespace Dumux {
template<class CouplingFacet, class Container, class IndexSet>
template<class CouplingFacet, class Container, class IndexSet, class Iterator>
class FacetIterator : public Dune::ForwardIteratorFacade<FacetIterator<CouplingFacet,
Container,
IndexSet>,
IndexSet,
Iterator>,
const CouplingFacet>
{
using ThisType = FacetIterator<CouplingFacet, Container, IndexSet>;
using Iterator = typename IndexSet::const_iterator;
using ThisType = FacetIterator<CouplingFacet, Container, IndexSet, Iterator>;
public:
FacetIterator(const Iterator& it, const Container& container)
: it_(it), container_(&container) {}
FacetIterator(const Iterator& it, const Container& container, const IndexSet& indices)
: it_(it), container_(&container), indices_(&indices) {}
FacetIterator() : it_(Iterator()), container_(nullptr) {}
FacetIterator() : it_(Iterator()), container_(nullptr), indices_(nullptr) {}
//! dereferencing yields a coupling facet
const CouplingFacet& dereference() const
......@@ -67,6 +67,21 @@ public:
return it_ == other.it_;
}
template<bool enable = std::is_integral<IndexSet>::value, std::enable_if_t<!enable, int> = 0>
void increment()
{
it_++;
if(it_ == (*indices_)[localIdx_].end())
{
localIdx_++;
auto it = std::find_if (indices_->begin()+localIdx_, indices_->end(), [] (const auto& idx) { return idx.size() > 0; });
if(it != indices_->end())
it_ = it->begin();
}
}
template<bool enable = std::is_integral<IndexSet>::value, std::enable_if_t<enable, int> = 0>
void increment()
{
it_++;
......@@ -75,8 +90,37 @@ public:
private:
Iterator it_;
const Container* container_;
const IndexSet* indices_;
std::size_t localIdx_{0};
};
template<class Mapper, std::size_t i>
auto couplingFacets(Dune::index_constant<i> domainI, const Mapper& mapper, std::size_t eIdx)
{
using FacetContainer = typename Mapper::FacetContainer;
using IndexContainerScvf = typename Mapper::IndexContainerScvf;
using IndexContainerElement = typename Mapper::IndexContainerElement;
using FacetIterator = FacetIterator<typename FacetContainer::value_type, FacetContainer, IndexContainerElement, typename IndexContainerScvf::const_iterator>;
const auto& indexSet = mapper.couplingFacetIdxMap(domainI).at(eIdx);
auto itBegin = std::find_if (indexSet.begin(), indexSet.end(), [] (const auto& idx) { return idx.size() > 0; });
auto itEnd = std::find_if(std::reverse_iterator(indexSet.end()),
std::reverse_iterator(indexSet.begin()),
[] (const auto& idx) { return idx.size() > 0; });
return Dune::IteratorRange<FacetIterator>(FacetIterator(itBegin->begin(), mapper.couplingFacets(), indexSet),
FacetIterator(itEnd->end(), mapper.couplingFacets(), indexSet));
}
template<class Mapper, std::size_t i>
auto couplingFacets(Dune::index_constant<i> domainI, const Mapper& mapper, std::size_t eIdx, std::size_t localScvfIdx)
{
using FacetContainer = typename Mapper::FacetContainer;
using IndexContainerScvf = typename Mapper::IndexContainerScvf;
using FacetIterator = FacetIterator<typename FacetContainer::value_type, FacetContainer, IndexContainerScvf, typename IndexContainerScvf::const_iterator>;
const auto& indexSet = mapper.couplingFacetIdxMap(domainI).at(eIdx)[localScvfIdx];
return Dune::IteratorRange<FacetIterator>(FacetIterator(indexSet->begin(), mapper.couplingFacets(), indexSet),
FacetIterator(indexSet->end(), mapper.couplingFacets(), indexSet));
}
/*!
* \ingroup StokesDarcyCoupling
......@@ -109,7 +153,6 @@ private:
static_assert(GridGeometry<porousMediumIdx>::discMethod == DiscretizationMethod::box, "The Darcy domain must use the Box discretization");
public:
// export the type describing a coupling segment
struct CouplingFacet
{
......@@ -124,14 +167,30 @@ public:
Geometry geometry;
};
template<class IndexSet>
inline Dune::IteratorRange<FacetIterator<CouplingFacet, std::vector<CouplingFacet>, IndexSet>>
couplingFacets(const IndexSet& indexSet)
{
using FacetIterator = FacetIterator<CouplingFacet, std::vector<CouplingFacet>, IndexSet>;
return Dune::IteratorRange<FacetIterator>(FacetIterator(indexSet.begin(), couplingFacets_),
FacetIterator(indexSet.end(), couplingFacets_));
}
using FacetContainer = std::vector<CouplingFacet>;
using IndexContainerScvf = std::vector<std::size_t>;
using IndexContainerElement = std::vector<IndexContainerScvf>;
// template<class IndexSet, bool enable = std::is_integral<IndexSet>::value, std::enable_if_t<enable, int> = 0>
// auto couplingFacets(const std::vector<IndexSet>& indexSet)
// {
// using FacetIterator = FacetIterator<CouplingFacet, FacetContainer, std::vector<IndexSet>, typename std::vector<IndexSet>::const_iterator>;
// return Dune::IteratorRange<FacetIterator>(FacetIterator(indexSet->begin(), couplingFacets_, indexSet),
// FacetIterator(indexSet->end(), couplingFacets_, indexSet));
// }
// template<class IndexSet, bool enable = std::is_integral<IndexSet>::value, std::enable_if_t<!enable, int> = 0>
// auto couplingFacets(const std::vector<IndexSet>& indexSet)
// {
// using FacetIterator = FacetIterator<CouplingFacet, FacetContainer, std::vector<IndexSet>, typename IndexSet::const_iterator>;
// auto itBegin = std::find_if (indexSet.begin(), indexSet.end(), [] (const auto& idx) { return idx.size() > 0; });
// auto itEnd = std::find_if(std::reverse_iterator(indexSet.end()),
// std::reverse_iterator(indexSet.begin()),
// [] (const auto& idx) { return idx.size() > 0; });
// return Dune::IteratorRange<FacetIterator>(FacetIterator(itBegin->begin(), couplingFacets_, indexSet),
// FacetIterator(itEnd->end(), couplingFacets_, indexSet));
// }
/*!
* \brief Main update routine
......@@ -151,48 +210,45 @@ public:
const auto& stokesFvGridGeometry = stokesProblem.gridGeometry();
const auto& darcyFvGridGeometry = darcyProblem.gridGeometry();
for(const auto& dataHandle : ffCouplingFacetIdxMap_)
for(const auto& elementData : ffCouplingFacetIdxMap_)
{
const auto stokesEIdx = dataHandle.first;
const auto stokesEIdx = elementData.first;
for (const auto& couplingData : dataHandle.second)
for(const auto& couplingFacet : Dumux::couplingFacets(Dune::index_constant<freeFlowIdx>(), *this, stokesEIdx))
{
for(const auto& couplingFacet : couplingFacets(couplingData))
{
const auto darcyEIdx = couplingFacet.pmEIdx;
const auto stokesScvfIdx = couplingFacet.ffScvfIdx;
const auto& stokesScvf = stokesFvGridGeometry.scvf(stokesScvfIdx);
const auto darcyEIdx = couplingFacet.pmEIdx;
const auto stokesScvfIdx = couplingFacet.ffScvfIdx;
const auto& stokesScvf = stokesFvGridGeometry.scvf(stokesScvfIdx);
const auto& darcyElement = darcyFvGridGeometry.element(darcyEIdx);
auto darcyFvGeometry = localView(darcyFvGridGeometry);
darcyFvGeometry.bind(darcyElement);
const auto& darcyElement = darcyFvGridGeometry.element(darcyEIdx);
auto darcyFvGeometry = localView(darcyFvGridGeometry);
darcyFvGeometry.bind(darcyElement);
const auto stokesElement = stokesFvGridGeometry.element(stokesEIdx);
auto stokesFvGeometry = localView(stokesFvGridGeometry);
stokesFvGeometry.bind(stokesElement);
const auto stokesElement = stokesFvGridGeometry.element(stokesEIdx);
auto stokesFvGeometry = localView(stokesFvGridGeometry);
stokesFvGeometry.bind(stokesElement);
darcyToStokesCellCenterStencils[darcyEIdx].push_back(stokesEIdx);
darcyToStokesFaceStencils[darcyEIdx].first.push_back(stokesScvf.dofIndex());
darcyToStokesFaceStencils[darcyEIdx].second.push_back(stokesScvf.index());
darcyToStokesCellCenterStencils[darcyEIdx].push_back(stokesEIdx);
darcyToStokesFaceStencils[darcyEIdx].first.push_back(stokesScvf.dofIndex());
darcyToStokesFaceStencils[darcyEIdx].second.push_back(stokesScvf.index());
for (auto&& scv : scvs(darcyFvGeometry))
{
stokesCellCenterToDarcyStencils[stokesEIdx].push_back(scv.dofIndex());
stokesFaceToDarcyStencils[stokesScvf.dofIndex()].push_back(scv.dofIndex());
}
for (auto&& scv : scvs(darcyFvGeometry))
{
stokesCellCenterToDarcyStencils[stokesEIdx].push_back(scv.dofIndex());
stokesFaceToDarcyStencils[stokesScvf.dofIndex()].push_back(scv.dofIndex());
}
if(!(slipCondition() == SlipCondition::BJS))
if(!(slipCondition() == SlipCondition::BJS))
{
const std::size_t numSubFaces = stokesScvf.pairData().size();
// Account for all interior sub-faces which include data from a boundary with slip condition
for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
{
const std::size_t numSubFaces = stokesScvf.pairData().size();
// Account for all interior sub-faces which include data from a boundary with slip condition
for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
{
const auto eIdx = stokesScvf.insideScvIdx();
const auto& lateralStokesScvf = stokesFvGeometry.scvf(eIdx, stokesScvf.pairData(localSubFaceIdx).localLateralFaceIdx);
if(lateralStokesScvf.dofIndex() != stokesScvf.dofIndex() && !lateralStokesScvf.boundary())
for (auto&& scv : scvs(darcyFvGeometry))
stokesFaceToDarcyStencils[lateralStokesScvf.dofIndex()].push_back(scv.dofIndex());
}
const auto eIdx = stokesScvf.insideScvIdx();
const auto& lateralStokesScvf = stokesFvGeometry.scvf(eIdx, stokesScvf.pairData(localSubFaceIdx).localLateralFaceIdx);
if(lateralStokesScvf.dofIndex() != stokesScvf.dofIndex() && !lateralStokesScvf.boundary())
for (auto&& scv : scvs(darcyFvGeometry))
stokesFaceToDarcyStencils[lateralStokesScvf.dofIndex()].push_back(scv.dofIndex());
}
}
}
......@@ -267,36 +323,37 @@ public:
*/
bool isCoupledDarcyScvf(std::size_t eIdx, std::size_t scvfLocalIdx) const
{
(pmCouplingFacetIdxMap_.count(eIdx) == 0) ? return false
: return pmCouplingFacetIdxMap_.at(eIdx)[scvfLocalIdx].size() > 0;
return (pmCouplingFacetIdxMap_.count(eIdx) == 0) ? false
: pmCouplingFacetIdxMap_.at(eIdx)[scvfLocalIdx].size() > 0;
}
/*!
* \brief A map that returns all Stokes elements coupled to a Darcy element
* \brief A map that returns all elements coupled to the other domain
*/
const auto& pmCouplingFacetIdxMap() const
template<std::size_t i>
const auto& couplingFacetIdxMap(Dune::index_constant<i> domainI) const
{
return pmCouplingFacetIdxMap_;
if constexpr (i == porousMediumIdx)
return pmCouplingFacetIdxMap_;
else if constexpr (i == freeFlowIdx)
return ffCouplingFacetIdxMap_;
}
/*!
* \brief A map that returns all Darcy elements coupled to a Stokes element
*/
const auto& ffCouplingFacetIdxMap() const
const CouplingFacet& couplingFacet(std::size_t idx) const
{
return ffCouplingFacetIdxMap_;
return couplingFacets_[idx];
}
const CouplingFacet& couplingFacet(std::size_t idx) const
const FacetContainer& couplingFacets() const
{
return couplingFacets_[idx];
return couplingFacets_;
}
private:
std::vector<CouplingFacet> couplingFacets_;
std::unordered_map<std::size_t, std::vector<std::vector<std::size_t>>> pmCouplingFacetIdxMap_;
std::unordered_map<std::size_t, std::vector<std::vector<std::size_t>>> ffCouplingFacetIdxMap_;
FacetContainer couplingFacets_;
std::unordered_map<std::size_t, IndexContainerElement> pmCouplingFacetIdxMap_;
std::unordered_map<std::size_t, IndexContainerElement> ffCouplingFacetIdxMap_;
};
} // end namespace Dumux
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment