24#ifndef DUMUX_CCMPFA_FACETCOUPLING_MANAGER_HH
25#define DUMUX_CCMPFA_FACETCOUPLING_MANAGER_HH
53template<
class MDTraits,
class CouplingMapper, std::
size_t bulkDomainId, std::
size_t lowDimDomainId>
55:
public FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainId, DiscretizationMethod::cctpfa>
60 using BulkIdType =
typename MDTraits::template SubDomain<bulkDomainId>::Index;
61 using LowDimIdType =
typename MDTraits::template SubDomain<lowDimDomainId>::Index;
62 static constexpr auto bulkId = BulkIdType();
63 static constexpr auto lowDimId = LowDimIdType();
66 template<std::
size_t id>
using SubDomainTypeTag =
typename MDTraits::template SubDomain<id>::TypeTag;
73 template<std::
size_t id>
using FVElementGeometry =
typename GridGeometry<id>::LocalView;
74 template<std::
size_t id>
using SubControlVolume =
typename GridGeometry<id>::SubControlVolume;
75 template<std::
size_t id>
using SubControlVolumeFace =
typename GridGeometry<id>::SubControlVolumeFace;
76 template<std::
size_t id>
using GridView =
typename GridGeometry<id>::GridView;
78 template<std::
size_t id>
using Element =
typename GridView<id>::template Codim<0>::Entity;
82 template<std::
size_t id>
using GridVolumeVariables =
typename GridVariables<id>::GridVolumeVariables;
83 template<std::
size_t id>
using ElementVolumeVariables =
typename GridVolumeVariables<id>::LocalView;
86 static constexpr int bulkDim = GridView<bulkDomainId>::dimension;
87 static constexpr int lowDimDim = GridView<lowDimDomainId>::dimension;
88 static constexpr auto bulkGridId = CouplingMapper::template gridId<bulkDim>();
89 static constexpr auto lowDimGridId = CouplingMapper::template gridId<lowDimDim>();
106 void init(std::shared_ptr< Problem<bulkId> > bulkProblem,
107 std::shared_ptr< Problem<lowDimId> > lowDimProblem,
108 std::shared_ptr< CouplingMapper > couplingMapper,
112 ParentType::init(bulkProblem, lowDimProblem, couplingMapper, curSol);
115 bulkScvfIsOnFacetElement_.assign(bulkProblem->gridGeometry().numScvf(),
false);
116 const auto& bulkMap = couplingMapper->couplingMap(bulkGridId, lowDimGridId);
117 for (
const auto& entry : bulkMap)
118 for (
const auto& couplingEntry : entry.second.elementToScvfMap)
119 for (
const auto& scvfIdx : couplingEntry.second)
120 bulkScvfIsOnFacetElement_[scvfIdx] =
true;
123 couplingMapperPtr_ = couplingMapper;
130 const SubControlVolumeFace<bulkId>& scvf)
const
131 {
return bulkScvfIsOnFacetElement_[scvf.index()]; }
140 template<
class TensorFunc>
142 const SubControlVolumeFace<bulkId>& scvf,
143 const TensorFunc& getT)
const
145 const auto lowDimElemIdx = this->getLowDimElementIndex(element, scvf);
147 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
148 const auto& s = map.find(this->bulkCouplingContext().elementIdx)->second.couplingElementStencil;
149 const auto& idxInContext = std::distance( s.begin(), std::find(s.begin(), s.end(), lowDimElemIdx) );
150 assert(std::find(s.begin(), s.end(), lowDimElemIdx) != s.end());
152 const auto& lowDimProblem = this->problem(lowDimId);
153 const auto& lowDimElement = lowDimProblem.gridGeometry().element(lowDimElemIdx);
154 const auto& lowDimVolVars = this->bulkCouplingContext().lowDimVolVars[idxInContext];
155 const auto& lowDimFvGeometry = this->bulkCouplingContext().lowDimFvGeometries[idxInContext];
158 return getT(lowDimProblem, lowDimElement, lowDimVolVars, lowDimFvGeometry, *scvs(lowDimFvGeometry).begin());
161 using ParentType::evalCouplingResidual;
168 template<
class LowDimLocalAssembler >
169 typename LocalResidual<lowDimId>::ElementResidualVector
171 const LowDimLocalAssembler& lowDimLocalAssembler,
173 GridIndexType<bulkId> dofIdxGlobalJ)
174 {
return evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId); }
181 template<
class LowDimLocalAssembler >
182 typename LocalResidual<lowDimId>::ElementResidualVector
186 assert(this->lowDimCouplingContext().isSet);
187 assert(this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element()) == this->lowDimCouplingContext().elementIdx);
190 typename LowDimLocalAssembler::LocalResidual::ElementResidualVector res(lowDimLocalAssembler.fvGeometry().numScv());
192 for (
const auto& scv : scvs(lowDimLocalAssembler.fvGeometry()))
193 res[scv.localDofIndex()] -= evalSourcesFromBulk(lowDimLocalAssembler.element(),
194 lowDimLocalAssembler.fvGeometry(),
195 lowDimLocalAssembler.curElemVolVars(),
204 const FVElementGeometry<lowDimId>& fvGeometry,
205 const ElementVolumeVariables<lowDimId>& elemVolVars,
206 const SubControlVolume<lowDimId>& scv)
209 assert(this->problem(lowDimId).gridGeometry().elementMapper().index(element) == this->lowDimCouplingContext().elementIdx);
211 NumEqVector<lowDimId> sources(0.0);
213 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
214 auto it = map.find(this->lowDimCouplingContext().elementIdx);
218 assert(this->lowDimCouplingContext().isSet);
219 for (
const auto& embedment : it->second.embedments)
224 const auto& coincidingScvfs = embedment.second;
225 const auto& scvfList = lowDimUsesBox ? std::vector<GridIndexType<lowDimId>>{ coincidingScvfs[scv.localDofIndex()] }
228 sources += this->evalBulkFluxes(this->problem(bulkId).gridGeometry().element(embedment.first),
229 *this->lowDimCouplingContext().bulkFvGeometry,
230 *this->lowDimCouplingContext().bulkElemVolVars,
231 *this->lowDimCouplingContext().bulkElemFluxVarsCache,
232 *this->lowDimCouplingContext().bulkLocalResidual,
243 template<
class JacobianPattern>
246 const auto& lowDimFVGridGeometry = this->problem(lowDimId).gridGeometry();
247 for (
const auto& element : elements(lowDimFVGridGeometry.gridView()))
250 const auto eIdx = lowDimFVGridGeometry.elementMapper().index(element);
251 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
252 auto it = map.find(eIdx);
258 const auto bulkElemIdx = it->second.embedments[0].first;
259 const auto& bulkMapEntry = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId).at(bulkElemIdx);
260 const auto& couplingStencil = bulkMapEntry.couplingStencil;
262 for (
auto globalJ : couplingStencil)
266 for (
int i = 0; i < element.subEntities(lowDimDim); ++i)
267 pattern.add(lowDimFVGridGeometry.vertexMapper().subIndex(element, i, lowDimDim), globalJ);
270 pattern.add(eIdx, globalJ);
277 template<
class JacobianPattern>
287 template<
class LowDimLocalAssembler,
class JacobianMatrixDiagBlock,
class Gr
idVariables>
289 const LowDimLocalAssembler& lowDimLocalAssembler,
290 const typename LowDimLocalAssembler::LocalResidual::ElementResidualVector&,
291 JacobianMatrixDiagBlock& A,
292 GridVariables& gridVariables)
296 if (!LowDimLocalAssembler::isImplicit())
300 auto updateContext = [&] (
auto elemIdx,
auto dofIdx,
auto priVars,
auto pvIdx)
303 auto& ldSol = this->curSol()[lowDimId];
304 ldSol[dofIdx][pvIdx] = priVars[pvIdx];
307 assert(this->bulkCouplingContext().isSet);
308 const auto& bulkMap = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
309 const auto& couplingElementStencil = bulkMap.find(this->bulkCouplingContext().elementIdx)->second.couplingElementStencil;
311 auto it = std::find(couplingElementStencil.begin(), couplingElementStencil.end(), elemIdx);
312 assert(it != couplingElementStencil.end());
313 const auto idxInContext = std::distance(couplingElementStencil.begin(), it);
315 auto& volVars = this->bulkCouplingContext().lowDimVolVars[idxInContext];
316 const auto& fvGeom = this->bulkCouplingContext().lowDimFvGeometries[idxInContext];
317 const auto& element = this->problem(lowDimId).gridGeometry().element(elemIdx);
322 const auto elemGeom = element.geometry();
327 volVars.update(
elementSolution(element, ldSol, this->problem(lowDimId).gridGeometry()),
328 this->problem(lowDimId),
330 fvGeom.scv(elemIdx) );
333 const auto contextElem = this->problem(bulkId).gridGeometry().element(this->bulkCouplingContext().elementIdx);
334 this->lowDimCouplingContext().bulkElemFluxVarsCache->update(contextElem,
335 *this->lowDimCouplingContext().bulkFvGeometry,
336 *this->lowDimCouplingContext().bulkElemVolVars);
339 const auto eIdx = this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element());
342 assert(this->lowDimCouplingContext().isSet);
343 assert(this->lowDimCouplingContext().elementIdx == eIdx);
346 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
347 auto it = map.find(eIdx);
349 evalLowDimSourceDerivatives_(updateContext, lowDimLocalAssembler, A);
353 template<
class LocalAssemblerI,
class JacobianMatrixDiagBlock,
class Gr
idVariables>
355 const LocalAssemblerI& localAssemblerI,
356 const typename LocalAssemblerI::LocalResidual::ElementResidualVector& origResiduals,
357 JacobianMatrixDiagBlock& A,
358 GridVariables& gridVariables)
363 template<
class UpdateContext,
class LowDimLocalAssembler,
class JacobianMatrixDiagBlock>
364 void evalLowDimSourceDerivatives_(
const UpdateContext& updateContext,
365 const LowDimLocalAssembler& lowDimLocalAssembler,
366 JacobianMatrixDiagBlock& A)
368 const auto& lowDimFVGridGeometry = this->problem(lowDimId).gridGeometry();
369 const auto eIdx = lowDimFVGridGeometry.elementMapper().index(lowDimLocalAssembler.element());
372 const auto bulkElemIdx = this->bulkCouplingContext().elementIdx;
373 const auto& bulkMapEntry = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId).at(bulkElemIdx);
374 const auto& couplingStencil = bulkMapEntry.couplingStencil;
375 const auto& couplingElementStencil = bulkMapEntry.couplingElementStencil;
378 const auto origResidual = evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId);
381 std::vector< std::decay_t<
decltype(couplingStencil[0])> > elemDofs;
382 elemDofs.reserve(lowDimLocalAssembler.fvGeometry().numScv());
383 for (
const auto& scv : scvs(lowDimLocalAssembler.fvGeometry()))
384 elemDofs.push_back(scv.dofIndex());
387 for (
const auto couplingElemIdx : couplingElementStencil)
390 if (couplingElemIdx == eIdx)
394 std::vector< std::decay_t<
decltype(couplingStencil[0])> > elemDofsJ;
397 const auto& elemJ = lowDimFVGridGeometry.element(couplingElemIdx);
398 for (
int i = 0; i < elemJ.subEntities(lowDimDim); ++i)
399 elemDofsJ.push_back(lowDimFVGridGeometry.vertexMapper().subIndex(elemJ, i, lowDimDim));
402 elemDofsJ.push_back(couplingElemIdx);
404 for (
auto dofIndex : elemDofsJ)
406 auto partialDerivs = origResidual;
407 const auto origPriVars = this->curSol()[lowDimId][dofIndex];
410 static constexpr auto numEq = std::decay_t<
decltype(origPriVars)>::dimension;
411 for (
int pvIdx = 0; pvIdx < numEq; pvIdx++)
416 auto evalResiduals = [&](Scalar<lowDimId> priVar)
418 auto priVars = origPriVars;
419 priVars[pvIdx] = priVar;
422 updateContext(couplingElemIdx, dofIndex, priVars, pvIdx);
423 return this->evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId);
426 static const int numDiffMethod = getParamFromGroup<int>(this->problem(lowDimId).paramGroup(),
"Assembly.NumericDifferenceMethod");
427 static const NumericEpsilon< Scalar<lowDimId>, numEq >
eps{this->problem(lowDimId).paramGroup()};
429 origResidual,
eps(origPriVars[pvIdx], pvIdx), numDiffMethod);
434 for (
const auto& scv : scvs(lowDimLocalAssembler.fvGeometry()))
435 for (
int eqIdx = 0; eqIdx < numEq; eqIdx++)
436 A[scv.dofIndex()][dofIndex][eqIdx][pvIdx] += partialDerivs[scv.indexInElement()][eqIdx];
439 updateContext(couplingElemIdx, dofIndex, origPriVars, pvIdx);
446 std::vector<bool> bulkScvfIsOnFacetElement_;
449 std::shared_ptr< CouplingMapper > couplingMapperPtr_;
An adapter class for local assemblers using numeric differentiation.
Defines the index types used for grid and local indices.
A class for numeric differentiation.
Element solution classes and factory functions.
The available discretization methods in Dumux.
DiscretizationMethod
The available discretization methods in Dumux.
Definition: method.hh:37
auto elementSolution(const Element &element, const SolutionVector &sol, const GridGeometry &gg) -> std::enable_if_t< GridGeometry::discMethod==DiscretizationMethod::box, BoxElementSolution< typename GridGeometry::LocalView, std::decay_t< decltype(std::declval< SolutionVector >()[0])> > >
Make an element solution for box schemes.
Definition: box/elementsolution.hh:115
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)
Free function that allows the creation of a volume variables object interpolated to a given position ...
Definition: multidomain/facet/couplingmanager.hh:52
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
constexpr double eps
epsilon for checking direction of scvf normals
Definition: test_tpfafvgeometry_nonconforming.cc:44
Struture to define the index types used for grid and local indices.
Definition: indextraits.hh:38
static void partialDerivative(const Function &function, Scalar x0, FunctionEvalType &derivative, const FunctionEvalType &fx0, const int numericDifferenceMethod=1)
Computes the derivative of a function with repect to a function parameter.
Definition: numericdifferentiation.hh:61
Property to specify the type of scalar values.
Definition: common/properties.hh:53
A vector of size number equations that can be used for Neumann fluxes, sources, residuals,...
Definition: common/properties.hh:61
Property to specify the type of a problem which has to be solved.
Definition: common/properties.hh:69
Definition: common/properties.hh:91
Definition: common/properties.hh:150
The grid variables object managing variable data on the grid (volvars/fluxvars cache)
Definition: common/properties.hh:190
void init(std::shared_ptr< Problem< bulkId > > bulkProblem, std::shared_ptr< Problem< lowDimId > > lowDimProblem, std::shared_ptr< CouplingMapper > couplingMapper, const SolutionVector &curSol)
Initialize the coupling manager.
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:106
bool isOnInteriorBoundary(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns true if a bulk scvf coincides with a facet element.
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:129
void extendJacobianPattern(LowDimIdType, JacobianPattern &pattern) const
Extend the jacobian pattern of the diagonal block of the lowdim domain by the elements that are in th...
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:244
typename MDTraits::SolutionVector SolutionVector
the type of the solution vector
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:96
LocalResidual< lowDimId >::ElementResidualVector evalCouplingResidual(LowDimIdType, const LowDimLocalAssembler &lowDimLocalAssembler, BulkIdType)
Evaluates the coupling element residual of a lower-dimensional domain element with respect to a dof i...
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:183
void evalAdditionalDomainDerivatives(LowDimIdType, const LowDimLocalAssembler &lowDimLocalAssembler, const typename LowDimLocalAssembler::LocalResidual::ElementResidualVector &, JacobianMatrixDiagBlock &A, GridVariables &gridVariables)
evaluate additional derivatives of the element residual of the low-dim domain with respect to dofs in...
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:288
void evalAdditionalDomainDerivatives(BulkIdType, const LocalAssemblerI &localAssemblerI, const typename LocalAssemblerI::LocalResidual::ElementResidualVector &origResiduals, JacobianMatrixDiagBlock &A, GridVariables &gridVariables)
The bulk domain has no additional derivatives.
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:354
LocalResidual< lowDimId >::ElementResidualVector evalCouplingResidual(LowDimIdType, const LowDimLocalAssembler &lowDimLocalAssembler, BulkIdType, GridIndexType< bulkId > dofIdxGlobalJ)
Evaluates the coupling element residual of a lower-dimensional domain element with respect to a dof i...
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:170
auto getLowDimTensor(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf, const TensorFunc &getT) const
returns the tensor within the low dim element coinciding with the given bulk scvf according to the pr...
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:141
void extendJacobianPattern(BulkIdType, JacobianPattern &pattern) const
The bulk domain has no extended jacobian pattern.
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:278
NumEqVector< lowDimId > evalSourcesFromBulk(const Element< lowDimId > &element, const FVElementGeometry< lowDimId > &fvGeometry, const ElementVolumeVariables< lowDimId > &elemVolVars, const SubControlVolume< lowDimId > &scv)
Computes the sources in a lower-dimensional sub-control volume stemming from the bulk domain.
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:203
Manages the coupling between bulk elements and lower dimensional elements where the coupling occurs a...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:53
Implementation for the coupling manager between two domains of dimension d and (d-1) for models consi...
Definition: multidomain/facet/couplingmanager.hh:95
Declares all properties used in Dumux.
The interface of the coupling manager for multi domain problems.