3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
multidomain/facet/box/couplingmanager.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_BOX_FACETCOUPLING_MANAGER_HH
25#define DUMUX_BOX_FACETCOUPLING_MANAGER_HH
26
27#include <algorithm>
28#include <cassert>
29
35
36namespace Dumux {
37
49template<class MDTraits, class CouplingMapper, std::size_t bulkDomainId, std::size_t lowDimDomainId>
50class FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainId, DiscretizationMethod::box>
51: public virtual CouplingManager< MDTraits >
52{
54
55 // convenience aliases and instances of the two domain ids
56 using BulkIdType = typename MDTraits::template SubDomain<bulkDomainId>::Index;
57 using LowDimIdType = typename MDTraits::template SubDomain<lowDimDomainId>::Index;
58 static constexpr auto bulkId = BulkIdType();
59 static constexpr auto lowDimId = LowDimIdType();
60
61 // the sub-domain type tags
62 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
63
64 // further types specific to the sub-problems
65 template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
66 template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
67 template<std::size_t id> using NumEqVector = GetPropType<SubDomainTypeTag<id>, Properties::NumEqVector>;
68 template<std::size_t id> using ElementBoundaryTypes = GetPropType<SubDomainTypeTag<id>, Properties::ElementBoundaryTypes>;
69 template<std::size_t id> using LocalResidual = GetPropType<SubDomainTypeTag<id>, Properties::LocalResidual>;
70
71 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
72 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
73 template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::SubControlVolume;
74 template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::SubControlVolumeFace;
75 template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
76 template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
77 template<std::size_t id> using GridIndexType = typename GridView<id>::IndexSet::IndexType;
78
79 template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
80 template<std::size_t id> using GridVolumeVariables = typename GridVariables<id>::GridVolumeVariables;
81 template<std::size_t id> using ElementVolumeVariables = typename GridVolumeVariables<id>::LocalView;
82 template<std::size_t id> using VolumeVariables = typename ElementVolumeVariables<id>::VolumeVariables;
83 template<std::size_t id> using GridFluxVariablesCache = typename GridVariables<id>::GridFluxVariablesCache;
84 template<std::size_t id> using ElementFluxVariablesCache = typename GridFluxVariablesCache<id>::LocalView;
85
86 // extract corresponding grid ids from the mapper
87 static constexpr int bulkDim = GridView<bulkDomainId>::dimension;
88 static constexpr int lowDimDim = GridView<lowDimDomainId>::dimension;
89 static constexpr auto bulkGridId = CouplingMapper::template gridId<bulkDim>();
90 static constexpr auto lowDimGridId = CouplingMapper::template gridId<lowDimDim>();
91
92 static constexpr bool lowDimUsesBox = GridGeometry<lowDimId>::discMethod == DiscretizationMethod::box;
93
100 struct BulkCouplingContext
101 {
102 bool isSet;
103 GridIndexType< bulkId > elementIdx;
104 std::vector< FVElementGeometry<lowDimId> > lowDimFvGeometries;
105 std::vector< std::vector<VolumeVariables<lowDimId>> > lowDimElemVolVars;
106
107 void reset()
108 {
109 lowDimFvGeometries.clear();
110 lowDimElemVolVars.clear();
111 isSet = false;
112 }
113 };
114
125 struct LowDimCouplingContext
126 {
127 private:
128 // For dim == dimworld in bulk domain, we know there is always going
129 // to be two bulk neighbors for a lower-dimensional element. Choose
130 // the container type depending on this
131 static constexpr int bulkDimWorld = GridView<bulkId>::dimensionworld;
132 static constexpr bool bulkIsSurfaceGrid = bulkDim != bulkDimWorld;
133
134 template< class T >
135 using ContainerType = typename std::conditional< bulkIsSurfaceGrid,
136 std::vector< T >,
137 std::array< T, 2 > >::type;
138
139 public:
140 bool isSet;
141 GridIndexType< lowDimId > elementIdx;
142 ContainerType< std::unique_ptr<FVElementGeometry<bulkId>> > bulkFvGeometries;
143 ContainerType< std::unique_ptr<ElementVolumeVariables<bulkId>> > bulkElemVolVars;
144 ContainerType< std::unique_ptr<ElementFluxVariablesCache<bulkId>> > bulkElemFluxVarsCache;
145 std::unique_ptr< LocalResidual<bulkId> > bulkLocalResidual;
146 ContainerType< ElementBoundaryTypes<bulkId> > bulkElemBcTypes;
147
148 void reset()
149 {
150 isSet = false;
151 auto resetPtr = [&] (auto&& uniquePtr) { uniquePtr.reset(); };
152 std::for_each(bulkFvGeometries.begin(), bulkFvGeometries.end(), resetPtr);
153 std::for_each(bulkElemVolVars.begin(), bulkElemVolVars.end(), resetPtr);
154 std::for_each(bulkElemFluxVarsCache.begin(), bulkElemFluxVarsCache.end(), resetPtr);
155 bulkLocalResidual.reset();
156 }
157
158 // resize containers for surface grids
159 template<bool s = bulkIsSurfaceGrid, std::enable_if_t<s, int> = 0>
160 void resize(std::size_t numEmbedments)
161 {
162 bulkFvGeometries.resize(numEmbedments);
163 bulkElemVolVars.resize(numEmbedments);
164 bulkElemFluxVarsCache.resize(numEmbedments);
165 bulkElemBcTypes.resize(numEmbedments);
166 }
167
168 // no memory reservation necessary for non-surface grids
169 template<bool s = bulkIsSurfaceGrid, std::enable_if_t<!s, int> = 0>
170 void resize(std::size_t numEmbedments)
171 {}
172 };
173
174public:
175
177 template<std::size_t i, std::size_t j = (i == bulkId) ? lowDimId : bulkId>
178 using CouplingStencilType = typename CouplingMapper::template Stencil< CouplingMapper::template gridId<GridView<j>::dimension>() >;
179
181 using SolutionVector = typename MDTraits::SolutionVector;
182
191 void init(std::shared_ptr< Problem<bulkId> > bulkProblem,
192 std::shared_ptr< Problem<lowDimId> > lowDimProblem,
193 std::shared_ptr< CouplingMapper > couplingMapper,
194 const SolutionVector& curSol)
195 {
196 couplingMapperPtr_ = couplingMapper;
197
198 // set the sub problems
199 this->setSubProblem(bulkProblem, bulkId);
200 this->setSubProblem(lowDimProblem, lowDimId);
201
202 // copy the solution vector
203 ParentType::updateSolution(curSol);
204
205 // determine all bulk elements that couple to low dim elements
206 const auto& bulkMap = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
207 bulkElemIsCoupled_.assign(bulkProblem->gridGeometry().gridView().size(0), false);
208 std::for_each( bulkMap.begin(),
209 bulkMap.end(),
210 [&] (const auto& entry) { bulkElemIsCoupled_[entry.first] = true; });
211 }
212
217 const Element<bulkId>& element,
218 LowDimIdType domainJ) const
219 {
220 const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
221
222 if (bulkElemIsCoupled_[eIdx])
223 {
224 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
225 auto it = map.find(eIdx);
226 assert(it != map.end());
227 return it->second.couplingStencil;
228 }
229
230 return getEmptyStencil(lowDimId);
231 }
232
237 const Element<lowDimId>& element,
238 BulkIdType domainJ) const
239 {
240 const auto eIdx = this->problem(lowDimId).gridGeometry().elementMapper().index(element);
241
242 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
243 auto it = map.find(eIdx);
244 if (it != map.end()) return it->second.couplingStencil;
245 else return getEmptyStencil(bulkId);
246 }
247
252 bool isCoupled(const Element<bulkId>& element,
253 const SubControlVolumeFace<bulkId>& scvf) const
254 { return scvf.interiorBoundary(); }
255
260 bool isOnInteriorBoundary(const Element<bulkId>& element,
261 const SubControlVolumeFace<bulkId>& scvf) const
262 { return isCoupled(element, scvf); }
263
267 const VolumeVariables<lowDimId>& getLowDimVolVars(const Element<bulkId>& element,
268 const SubControlVolumeFace<bulkId>& scvf) const
269 {
270 const auto eIdx = this->problem(bulkId).gridGeometry().elementMapper().index(element);
271 assert(bulkContext_.isSet);
272 assert(bulkElemIsCoupled_[eIdx]);
273
274 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
275 const auto& couplingData = map.find(eIdx)->second;
276
277 // search the low dim element idx this scvf is embedded in
278 // and find out the local index of the scv it couples to
279 unsigned int coupledScvIdx;
280 auto it = std::find_if( couplingData.elementToScvfMap.begin(),
281 couplingData.elementToScvfMap.end(),
282 [&] (auto& dataPair)
283 {
284 const auto& scvfList = dataPair.second;
285 auto it = std::find(scvfList.begin(), scvfList.end(), scvf.index());
286 coupledScvIdx = std::distance(scvfList.begin(), it);
287 return it != scvfList.end();
288 } );
289
290 assert(it != couplingData.elementToScvfMap.end());
291 const auto lowDimElemIdx = it->first;
292 const auto& s = map.find(bulkContext_.elementIdx)->second.couplingElementStencil;
293 const auto& idxInContext = std::distance( s.begin(), std::find(s.begin(), s.end(), lowDimElemIdx) );
294 assert(std::find(s.begin(), s.end(), lowDimElemIdx) != s.end());
295
296 if (lowDimUsesBox)
297 return bulkContext_.lowDimElemVolVars[idxInContext][coupledScvIdx];
298 else
299 return bulkContext_.lowDimElemVolVars[idxInContext][0];
300 }
301
305 const Element<lowDimId> getLowDimElement(const Element<bulkId>& element,
306 const SubControlVolumeFace<bulkId>& scvf) const
307 {
308 const auto lowDimElemIdx = getLowDimElementIndex(element, scvf);
309 return this->problem(lowDimId).gridGeometry().element(lowDimElemIdx);
310 }
311
315 const GridIndexType<lowDimId> getLowDimElementIndex(const Element<bulkId>& element,
316 const SubControlVolumeFace<bulkId>& scvf) const
317 {
318 const auto eIdx = this->problem(bulkId).gridGeometry().elementMapper().index(element);
319
320 assert(bulkElemIsCoupled_[eIdx]);
321 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
322 const auto& couplingData = map.at(eIdx);
323
324 // search the low dim element idx this scvf is embedded in
325 auto it = std::find_if( couplingData.elementToScvfMap.begin(),
326 couplingData.elementToScvfMap.end(),
327 [&] (auto& dataPair)
328 {
329 const auto& scvfList = dataPair.second;
330 auto it = std::find(scvfList.begin(), scvfList.end(), scvf.index());
331 return it != scvfList.end();
332 } );
333
334 assert(it != couplingData.elementToScvfMap.end());
335 return it->first;
336 }
337
344 template< class BulkLocalAssembler >
345 typename LocalResidual<bulkId>::ElementResidualVector
347 const BulkLocalAssembler& bulkLocalAssembler,
348 LowDimIdType,
349 GridIndexType<lowDimId> dofIdxGlobalJ)
350 {
351 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
352
353 assert(bulkContext_.isSet);
354 assert(bulkElemIsCoupled_[bulkContext_.elementIdx]);
355 assert(map.find(bulkContext_.elementIdx) != map.end());
356 assert(bulkContext_.elementIdx == this->problem(bulkId).gridGeometry().elementMapper().index(bulkLocalAssembler.element()));
357
358 const auto& fvGeometry = bulkLocalAssembler.fvGeometry();
359 typename LocalResidual<bulkId>::ElementResidualVector res(fvGeometry.numScv());
360 res = 0.0;
361
362 // compute fluxes across the coupling scvfs
363 const auto& couplingScvfs = map.find(bulkContext_.elementIdx)->second.dofToCouplingScvfMap.at(dofIdxGlobalJ);
364 for (auto scvfIdx : couplingScvfs)
365 {
366 const auto& scvf = fvGeometry.scvf(scvfIdx);
367 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
368 res[insideScv.localDofIndex()] += evalBulkFluxes_( bulkLocalAssembler.element(),
369 bulkLocalAssembler.fvGeometry(),
370 bulkLocalAssembler.curElemVolVars(),
371 bulkLocalAssembler.elemBcTypes(),
372 bulkLocalAssembler.elemFluxVarsCache(),
373 bulkLocalAssembler.localResidual(),
374 std::array<GridIndexType<bulkId>, 1>({scvfIdx}) );
375 }
376
377 return res;
378 }
379
386 template< class LowDimLocalAssembler >
387 typename LocalResidual<lowDimId>::ElementResidualVector
389 const LowDimLocalAssembler& lowDimLocalAssembler,
390 BulkIdType,
391 GridIndexType<bulkId> dofIdxGlobalJ)
392 {
393 // make sure this is called for the element for which the context was set
394 assert(lowDimContext_.isSet);
395 assert(this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element()) == lowDimContext_.elementIdx);
396
397 // evaluate sources for all scvs in lower-dimensional element
398 typename LocalResidual<lowDimId>::ElementResidualVector res(lowDimLocalAssembler.fvGeometry().numScv());
399 res = 0.0;
400 for (const auto& scv : scvs(lowDimLocalAssembler.fvGeometry()))
401 res[scv.localDofIndex()] -= evalSourcesFromBulk(lowDimLocalAssembler.element(),
402 lowDimLocalAssembler.fvGeometry(),
403 lowDimLocalAssembler.curElemVolVars(),
404 scv);
405 return res;
406 }
407
411 NumEqVector<lowDimId> evalSourcesFromBulk(const Element<lowDimId>& element,
412 const FVElementGeometry<lowDimId>& fvGeometry,
413 const ElementVolumeVariables<lowDimId>& elemVolVars,
414 const SubControlVolume<lowDimId>& scv)
415 {
416 // make sure the this is called for the element of the context
417 assert(this->problem(lowDimId).gridGeometry().elementMapper().index(element) == lowDimContext_.elementIdx);
418
419 NumEqVector<lowDimId> sources(0.0);
420 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
421 auto it = map.find(lowDimContext_.elementIdx);
422 if (it == map.end())
423 return sources;
424
425 assert(lowDimContext_.isSet);
426 for (unsigned int i = 0; i < it->second.embedments.size(); ++i)
427 {
428 const auto& embedment = it->second.embedments[i];
429
430 // list of scvfs in the bulk domain whose fluxes enter this scv
431 // if low dim domain uses tpfa, this is all scvfs lying on this element
432 // if it uses box, it is the one scvf coinciding with the given scv
433 const auto& coincidingScvfs = embedment.second;
434 const auto& scvfList = lowDimUsesBox ? std::vector<GridIndexType<lowDimId>>{ coincidingScvfs[scv.localDofIndex()] }
435 : coincidingScvfs;
436
437 sources += evalBulkFluxes_(this->problem(bulkId).gridGeometry().element(embedment.first),
438 *(lowDimContext_.bulkFvGeometries[i]),
439 *(lowDimContext_.bulkElemVolVars[i]),
440 lowDimContext_.bulkElemBcTypes[i],
441 *(lowDimContext_.bulkElemFluxVarsCache[i]),
442 *lowDimContext_.bulkLocalResidual,
443 scvfList);
444 }
445
446 return sources;
447 }
448
454 template< class Assembler >
455 void bindCouplingContext(BulkIdType, const Element<bulkId>& element, const Assembler& assembler)
456 {
457 // clear context
458 bulkContext_.reset();
459
460 // set index in context in any case
461 const auto bulkElemIdx = this->problem(bulkId).gridGeometry().elementMapper().index(element);
462 bulkContext_.elementIdx = bulkElemIdx;
463
464 // if element is coupled, actually set the context
465 if (bulkElemIsCoupled_[bulkElemIdx])
466 {
467 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
468
469 auto it = map.find(bulkElemIdx); assert(it != map.end());
470 const auto& elementStencil = it->second.couplingElementStencil;
471 bulkContext_.lowDimFvGeometries.reserve(elementStencil.size());
472 bulkContext_.lowDimElemVolVars.reserve(elementStencil.size());
473
474 // local view on the bulk fv geometry
475 auto bulkFvGeometry = localView(this->problem(bulkId).gridGeometry());
476 bulkFvGeometry.bindElement(element);
477
478 for (const auto lowDimElemIdx : elementStencil)
479 {
480 const auto& ldGridGeometry = this->problem(lowDimId).gridGeometry();
481
482 const auto& ldSol = Assembler::isImplicit() ? this->curSol()[lowDimId] : assembler.prevSol()[lowDimId];
483 const auto elemJ = ldGridGeometry.element(lowDimElemIdx);
484 auto fvGeom = localView(ldGridGeometry);
485 fvGeom.bindElement(elemJ);
486
487 std::vector<VolumeVariables<lowDimId>> elemVolVars(fvGeom.numScv());
488 const auto& coupledScvfIndices = it->second.elementToScvfMap.at(lowDimElemIdx);
489 makeCoupledLowDimElemVolVars_(element, bulkFvGeometry, elemJ, fvGeom, ldSol, coupledScvfIndices, elemVolVars);
490
491 bulkContext_.isSet = true;
492 bulkContext_.lowDimFvGeometries.emplace_back( std::move(fvGeom) );
493 bulkContext_.lowDimElemVolVars.emplace_back( std::move(elemVolVars) );
494 }
495 }
496 }
497
507 template< class Assembler >
508 void bindCouplingContext(LowDimIdType, const Element<lowDimId>& element, const Assembler& assembler)
509 {
510 // reset contexts
511 bulkContext_.reset();
512 lowDimContext_.reset();
513
514 // set index in context in any case
515 const auto lowDimElemIdx = this->problem(lowDimId).gridGeometry().elementMapper().index(element);
516 lowDimContext_.elementIdx = lowDimElemIdx;
517
518 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
519 auto it = map.find(lowDimElemIdx);
520
521 // if element is coupled, actually set the context
522 if (it != map.end())
523 {
524 // first bind the low dim context for the first neighboring bulk element
525 // this will set up the lower-dimensional data on this current lowdim element
526 // which will be enough to compute the fluxes in all neighboring elements
527 const auto& bulkGridGeom = this->problem(bulkId).gridGeometry();
528 const auto bulkElem = bulkGridGeom.element(it->second.embedments[0].first);
529 bindCouplingContext(bulkId, bulkElem, assembler);
530
531 // reserve space in the context and bind local views of neighbors
532 const auto& embedments = it->second.embedments;
533 const auto numEmbedments = embedments.size();
534 lowDimContext_.resize(numEmbedments);
535 for (unsigned int i = 0; i < numEmbedments; ++i)
536 {
537 auto bulkFvGeom = localView(bulkGridGeom);
538 auto bulkElemVolVars = localView(assembler.gridVariables(bulkId).curGridVolVars());
539 auto bulkElemFluxVarsCache = localView(assembler.gridVariables(bulkId).gridFluxVarsCache());
540
541 const auto& bulkSol = Assembler::isImplicit() ? this->curSol()[bulkId] : assembler.prevSol()[bulkId];
542 const auto curBulkElem = bulkGridGeom.element(embedments[i].first);
543 bulkFvGeom.bind(curBulkElem);
544 bulkElemVolVars.bind(curBulkElem, bulkFvGeom, bulkSol);
545 bulkElemFluxVarsCache.bind(curBulkElem, bulkFvGeom, bulkElemVolVars);
546
547 lowDimContext_.isSet = true;
548 lowDimContext_.bulkElemBcTypes[i].update(this->problem(bulkId), curBulkElem, bulkFvGeom);
549 lowDimContext_.bulkFvGeometries[i] = std::make_unique< FVElementGeometry<bulkId> >( std::move(bulkFvGeom) );
550 lowDimContext_.bulkElemVolVars[i] = std::make_unique< ElementVolumeVariables<bulkId> >( std::move(bulkElemVolVars) );
551 lowDimContext_.bulkElemFluxVarsCache[i] = std::make_unique< ElementFluxVariablesCache<bulkId> >( std::move(bulkElemFluxVarsCache) );
552 }
553
554 // finally, set the local residual
555 lowDimContext_.bulkLocalResidual = std::make_unique< LocalResidual<bulkId> >(assembler.localResidual(bulkId));
556 }
557 }
558
563 template< class BulkLocalAssembler >
564 void updateCouplingContext(BulkIdType domainI,
565 const BulkLocalAssembler& bulkLocalAssembler,
566 LowDimIdType domainJ,
567 GridIndexType<lowDimId> dofIdxGlobalJ,
568 const PrimaryVariables<lowDimId>& priVarsJ,
569 unsigned int pvIdxJ)
570 {
571 // communicate deflected solution
572 ParentType::updateCouplingContext(domainI, bulkLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
573
574 // Since coupling only occurs via the fluxes, the context does not
575 // have to be updated in explicit time discretization schemes, where
576 // they are strictly evaluated on the old time level
577 if (!BulkLocalAssembler::isImplicit())
578 return;
579
580 // skip the rest if context is empty
581 if (bulkContext_.isSet)
582 {
583 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
584 const auto& couplingEntry = map.at(bulkContext_.elementIdx);
585 const auto& couplingElemStencil = couplingEntry.couplingElementStencil;
586 const auto& ldGridGeometry = this->problem(lowDimId).gridGeometry();
587
588 // find the low-dim elements in coupling stencil, where this dof is contained in
589 const auto couplingElements = [&] ()
590 {
591 if (lowDimUsesBox)
592 {
593 std::vector< Element<lowDimId> > lowDimElems;
594 std::for_each( couplingElemStencil.begin(), couplingElemStencil.end(),
595 [&] (auto lowDimElemIdx)
596 {
597 auto element = ldGridGeometry.element(lowDimElemIdx);
598 for (unsigned int i = 0; i < element.geometry().corners(); ++i)
599 {
600 const auto dofIdx = ldGridGeometry.vertexMapper().subIndex(element, i, lowDimDim);
601 if (dofIdxGlobalJ == dofIdx) { lowDimElems.emplace_back( std::move(element) ); break; }
602 }
603 } );
604 return lowDimElems;
605 }
606 // dof index = element index for cc schemes
607 else
608 return std::vector<Element<lowDimId>>( {ldGridGeometry.element(dofIdxGlobalJ)} );
609 } ();
610
611 // update all necessary vol vars in context
612 for (const auto& element : couplingElements)
613 {
614 // find index in coupling context
615 const auto eIdxGlobal = ldGridGeometry.elementMapper().index(element);
616 auto it = std::find(couplingElemStencil.begin(), couplingElemStencil.end(), eIdxGlobal);
617 const auto idxInContext = std::distance(couplingElemStencil.begin(), it);
618 assert(it != couplingElemStencil.end());
619
620 auto& elemVolVars = bulkContext_.lowDimElemVolVars[idxInContext];
621 const auto& fvGeom = bulkContext_.lowDimFvGeometries[idxInContext];
622 const auto& coupledScvfIndices = couplingEntry.elementToScvfMap.at(eIdxGlobal);
623 makeCoupledLowDimElemVolVars_(bulkLocalAssembler.element(), bulkLocalAssembler.fvGeometry(),
624 element, fvGeom, this->curSol()[lowDimId], coupledScvfIndices, elemVolVars);
625 }
626 }
627 }
628
633 template< class BulkLocalAssembler >
634 void updateCouplingContext(BulkIdType domainI,
635 const BulkLocalAssembler& bulkLocalAssembler,
636 BulkIdType domainJ,
637 GridIndexType<bulkId> dofIdxGlobalJ,
638 const PrimaryVariables<bulkId>& priVarsJ,
639 unsigned int pvIdxJ)
640 {
641 // communicate deflected solution
642 ParentType::updateCouplingContext(domainI, bulkLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
643 }
644
650 template< class LowDimLocalAssembler >
651 void updateCouplingContext(LowDimIdType domainI,
652 const LowDimLocalAssembler& lowDimLocalAssembler,
653 BulkIdType domainJ,
654 GridIndexType<bulkId> dofIdxGlobalJ,
655 const PrimaryVariables<bulkId>& priVarsJ,
656 unsigned int pvIdxJ)
657 {
658 // communicate deflected solution
659 ParentType::updateCouplingContext(domainI, lowDimLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
660
661 // Since coupling only occurs via the fluxes, the context does not
662 // have to be updated in explicit time discretization schemes, where
663 // they are strictly evaluated on the old time level
664 if (!LowDimLocalAssembler::isImplicit())
665 return;
666
667 // skip the rest if context is empty
668 if (lowDimContext_.isSet)
669 {
670 assert(lowDimContext_.elementIdx == this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element()));
671
672 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
673 auto it = map.find(lowDimContext_.elementIdx);
674
675 assert(it != map.end());
676 const auto& embedments = it->second.embedments;
677 const auto& bulkGridGeometry = this->problem(bulkId).gridGeometry();
678
679 // TODO more efficient (only one dof update per embedment)
680 // update the elem volvars in context of those elements that contain globalJ
681 unsigned int embedmentIdx = 0;
682 for (const auto& embedment : embedments)
683 {
684 const auto elementJ = bulkGridGeometry.element(embedment.first);
685
686 for (unsigned int i = 0; i < elementJ.subEntities(bulkDim); ++i)
687 {
688 const auto dofIdx = bulkGridGeometry.vertexMapper().subIndex(elementJ, i, bulkDim);
689 if (dofIdx == dofIdxGlobalJ)
690 {
691 // element contains the deflected dof
692 const auto& fvGeom = *lowDimContext_.bulkFvGeometries[embedmentIdx];
693 (*lowDimContext_.bulkElemVolVars[embedmentIdx]).bindElement(elementJ, fvGeom, this->curSol()[bulkId]);
694 }
695 }
696
697 // keep track of embedments
698 embedmentIdx++;
699 }
700 }
701 }
702
709 template< class LowDimLocalAssembler >
710 void updateCouplingContext(LowDimIdType domainI,
711 const LowDimLocalAssembler& lowDimLocalAssembler,
712 LowDimIdType domainJ,
713 GridIndexType<lowDimId> dofIdxGlobalJ,
714 const PrimaryVariables<lowDimId>& priVarsJ,
715 unsigned int pvIdxJ)
716 {
717 // communicate deflected solution
718 ParentType::updateCouplingContext(domainI, lowDimLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
719
720 // Since coupling only occurs via the fluxes, the context does not
721 // have to be updated in explicit time discretization schemes, where
722 // they are strictly evaluated on the old time level
723 if (!LowDimLocalAssembler::isImplicit())
724 return;
725
726 // skip the rest if context is empty
727 if (lowDimContext_.isSet)
728 {
729 assert(bulkContext_.isSet);
730 assert(lowDimContext_.elementIdx == this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element()));
731
732 // update the corresponding vol vars in the bulk context
733 const auto& bulkMap = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
734 const auto& bulkCouplingEntry = bulkMap.at(bulkContext_.elementIdx);
735 const auto& couplingElementStencil = bulkCouplingEntry.couplingElementStencil;
736 auto it = std::find(couplingElementStencil.begin(), couplingElementStencil.end(), lowDimContext_.elementIdx);
737 assert(it != couplingElementStencil.end());
738 const auto idxInContext = std::distance(couplingElementStencil.begin(), it);
739
740 // get neighboring bulk element from the bulk context (is the same elemet as first entry in low dim context)
741 const auto& bulkElement = this->problem(bulkId).gridGeometry().element(bulkContext_.elementIdx);
742 const auto& bulkFvGeometry = *lowDimContext_.bulkFvGeometries[0];
743
744 auto& elemVolVars = bulkContext_.lowDimElemVolVars[idxInContext];
745 const auto& fvGeom = bulkContext_.lowDimFvGeometries[idxInContext];
746 const auto& element = lowDimLocalAssembler.element();
747 const auto& scvfIndices = bulkCouplingEntry.elementToScvfMap.at(lowDimContext_.elementIdx);
748
749 makeCoupledLowDimElemVolVars_(bulkElement, bulkFvGeometry, element, fvGeom,
750 this->curSol()[lowDimId], scvfIndices, elemVolVars);
751 }
752 }
753
755 template<std::size_t id, std::enable_if_t<(id == bulkId || id == lowDimId), int> = 0>
756 const typename CouplingMapper::template Stencil<id>&
757 getEmptyStencil(Dune::index_constant<id>) const
758 { return std::get<(id == bulkId ? 0 : 1)>(emptyStencilTuple_); }
759
760private:
761
766 template<class LowDimSolution, class ScvfIdxStorage>
767 void makeCoupledLowDimElemVolVars_(const Element<bulkId>& bulkElement,
768 const FVElementGeometry<bulkId>& bulkFvGeometry,
769 const Element<lowDimId>& lowDimElement,
770 const FVElementGeometry<lowDimId>& lowDimFvGeometry,
771 const LowDimSolution& lowDimSol,
772 const ScvfIdxStorage& coupledBulkScvfIndices,
773 std::vector<VolumeVariables<lowDimId>>& elemVolVars)
774 {
775 const auto lowDimElemSol = elementSolution(lowDimElement, lowDimSol, lowDimFvGeometry.gridGeometry());
776
777 // for tpfa, make only one volume variables object for the element
778 // for the update, use the first (and only - for tpfa) scv of the element
779 if (!lowDimUsesBox)
780 elemVolVars[0].update(lowDimElemSol, this->problem(lowDimId), lowDimElement, *scvs(lowDimFvGeometry).begin());
781 // for box, make vol vars either:
782 // - in each scv center, which geometrically coincides with the scvf integration point (Neumann coupling)
783 // - at the vertex position (Dirichlet coupling)
784 else
785 {
786 // recall that the scvfs in the coupling map were ordered such
787 // that the i-th scvf in the map couples to the i-th lowdim scv
788 for (unsigned int i = 0; i < coupledBulkScvfIndices.size(); ++i)
789 {
790 const auto& scvf = bulkFvGeometry.scvf(coupledBulkScvfIndices[i]);
791 const auto bcTypes = this->problem(bulkId).interiorBoundaryTypes(bulkElement, scvf);
792 if (bcTypes.hasOnlyNeumann())
793 FacetCoupling::makeInterpolatedVolVars(elemVolVars[i], this->problem(lowDimId), lowDimSol,
794 lowDimFvGeometry, lowDimElement, lowDimElement.geometry(),
795 scvf.ipGlobal());
796 else
797 elemVolVars[i].update( lowDimElemSol, this->problem(lowDimId), lowDimElement, lowDimFvGeometry.scv(i) );
798 }
799 }
800 }
801
803 template<class BulkScvfIndices>
804 NumEqVector<bulkId> evalBulkFluxes_(const Element<bulkId>& elementI,
805 const FVElementGeometry<bulkId>& fvGeometry,
806 const ElementVolumeVariables<bulkId>& elemVolVars,
807 const ElementBoundaryTypes<bulkId>& elemBcTypes,
808 const ElementFluxVariablesCache<bulkId>& elemFluxVarsCache,
809 const LocalResidual<bulkId>& localResidual,
810 const BulkScvfIndices& scvfIndices) const
811 {
812
813 NumEqVector<bulkId> coupledFluxes(0.0);
814 for (const auto& scvfIdx : scvfIndices)
815 coupledFluxes += localResidual.evalFlux(this->problem(bulkId),
816 elementI,
817 fvGeometry,
818 elemVolVars,
819 elemBcTypes,
820 elemFluxVarsCache,
821 fvGeometry.scvf(scvfIdx));
822 return coupledFluxes;
823 }
824
825 std::shared_ptr<CouplingMapper> couplingMapperPtr_;
826
829 std::vector<bool> bulkElemIsCoupled_;
830
832 using BulkStencil = typename CouplingMapper::template Stencil<bulkId>;
833 using LowDimStencil = typename CouplingMapper::template Stencil<lowDimId>;
834 std::tuple<BulkStencil, LowDimStencil> emptyStencilTuple_;
835
837 BulkCouplingContext bulkContext_;
838 LowDimCouplingContext lowDimContext_;
839};
840
841} // end namespace Dumux
842
843#endif
Element solution classes and factory functions.
The available discretization methods in Dumux.
ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: geometry/distance.hh:138
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
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
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
A vector of primary variables.
Definition: common/properties.hh:49
A vector of size number equations that can be used for Neumann fluxes, sources, residuals,...
Definition: common/properties.hh:51
Property to specify the type of a problem which has to be solved.
Definition: common/properties.hh:57
Definition: common/properties.hh:77
Stores the boundary types on an element.
Definition: common/properties.hh:98
Definition: common/properties.hh:101
The grid variables object managing variable data on the grid (volvars/fluxvars cache)
Definition: common/properties.hh:122
Definition: multidomain/couplingmanager.hh:46
void updateCouplingContext(LowDimIdType domainI, const LowDimLocalAssembler &lowDimLocalAssembler, LowDimIdType domainJ, GridIndexType< lowDimId > dofIdxGlobalJ, const PrimaryVariables< lowDimId > &priVarsJ, unsigned int pvIdxJ)
After deflecting the solution of the lower-dimensional domain has been deflected during the assembly ...
Definition: multidomain/facet/box/couplingmanager.hh:710
const GridIndexType< lowDimId > getLowDimElementIndex(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns the index of the lower-dimensional element coinciding with a bulk scvf.
Definition: multidomain/facet/box/couplingmanager.hh:315
bool isCoupled(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns true if a bulk scvf flux depends on data in the facet domain.
Definition: multidomain/facet/box/couplingmanager.hh:252
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 element stemming from the bulk domain.
Definition: multidomain/facet/box/couplingmanager.hh:411
LocalResidual< bulkId >::ElementResidualVector evalCouplingResidual(BulkIdType, const BulkLocalAssembler &bulkLocalAssembler, LowDimIdType, GridIndexType< lowDimId > dofIdxGlobalJ)
Evaluates the coupling element residual of a bulk domain element with respect to a dof in the lower-d...
Definition: multidomain/facet/box/couplingmanager.hh:346
const CouplingStencilType< lowDimId > & couplingStencil(LowDimIdType domainI, const Element< lowDimId > &element, BulkIdType domainJ) const
The coupling stencil of the lower-dimensional domain with the bulk domain.
Definition: multidomain/facet/box/couplingmanager.hh:236
void bindCouplingContext(BulkIdType, const Element< bulkId > &element, const Assembler &assembler)
For the assembly of the element residual of a bulk domain element we need to prepare all variables of...
Definition: multidomain/facet/box/couplingmanager.hh:455
void updateCouplingContext(LowDimIdType domainI, const LowDimLocalAssembler &lowDimLocalAssembler, BulkIdType domainJ, GridIndexType< bulkId > dofIdxGlobalJ, const PrimaryVariables< bulkId > &priVarsJ, unsigned int pvIdxJ)
After deflecting the solution of the bulk domain, we have to update the element volume variables of t...
Definition: multidomain/facet/box/couplingmanager.hh:651
void updateCouplingContext(BulkIdType domainI, const BulkLocalAssembler &bulkLocalAssembler, BulkIdType domainJ, GridIndexType< bulkId > dofIdxGlobalJ, const PrimaryVariables< bulkId > &priVarsJ, unsigned int pvIdxJ)
Update the coupling context for a derivative bulk -> bulk. Here, we simply have to update the solutio...
Definition: multidomain/facet/box/couplingmanager.hh:634
void updateCouplingContext(BulkIdType domainI, const BulkLocalAssembler &bulkLocalAssembler, LowDimIdType domainJ, GridIndexType< lowDimId > dofIdxGlobalJ, const PrimaryVariables< lowDimId > &priVarsJ, unsigned int pvIdxJ)
After deflecting the solution of the lower-dimensional domain, we have to update the element volume v...
Definition: multidomain/facet/box/couplingmanager.hh:564
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/box/couplingmanager.hh:388
const Element< lowDimId > getLowDimElement(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns the lower-dimensional element coinciding with a bulk scvf.
Definition: multidomain/facet/box/couplingmanager.hh:305
const VolumeVariables< lowDimId > & getLowDimVolVars(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns the vol vars of a lower-dimensional element coinciding with a bulk scvf.
Definition: multidomain/facet/box/couplingmanager.hh:267
const CouplingStencilType< bulkId > & couplingStencil(BulkIdType domainI, const Element< bulkId > &element, LowDimIdType domainJ) const
The coupling stencil of a given bulk domain element.
Definition: multidomain/facet/box/couplingmanager.hh:216
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/box/couplingmanager.hh:260
typename CouplingMapper::template Stencil< CouplingMapper::template gridId< GridView< j >::dimension >() > CouplingStencilType
types used for coupling stencils
Definition: multidomain/facet/box/couplingmanager.hh:178
const CouplingMapper::template Stencil< id > & getEmptyStencil(Dune::index_constant< id >) const
Empty stencil to be returned for elements that aren't coupled.
Definition: multidomain/facet/box/couplingmanager.hh:757
void bindCouplingContext(LowDimIdType, const Element< lowDimId > &element, const Assembler &assembler)
For the assembly of the element residual of a bulk domain element we need to prepare the local views ...
Definition: multidomain/facet/box/couplingmanager.hh:508
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/box/couplingmanager.hh:191
typename MDTraits::SolutionVector SolutionVector
the type of the solution vector
Definition: multidomain/facet/box/couplingmanager.hh:181
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.