couplingmanager.hh 27.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
 * \ingroup StokesDarcyCoupling
 * \copydoc Dumux::StokesDarcyCouplingManager
 */

#ifndef DUMUX_STOKES_DARCY_BOX_COUPLINGMANAGER_HH
#define DUMUX_STOKES_DARCY_BOX_COUPLINGMANAGER_HH

#include <utility>
#include <memory>

#include <dune/common/float_cmp.hh>
#include <dune/common/exceptions.hh>
#include <dumux/common/properties.hh>
#include <dumux/multidomain/staggeredcouplingmanager.hh>
#include <dumux/discretization/staggered/elementsolution.hh>

37
38
#include "../couplingmanager.hh"
#include "../couplingdata.hh"
39
40
41
42
43
44
45
46
47
48
#include "couplingmapper.hh"

namespace Dumux {

/*!
 * \ingroup StokesDarcyCoupling
 * \brief Coupling manager for Stokes and Darcy domains with equal dimension.
 */
template<class MDTraits>
class StokesDarcyCouplingManagerImplementation<MDTraits, DiscretizationMethod::box>
49
: public virtual StaggeredCouplingManager<MDTraits>, public FreeFlowPorousMediumCouplingManagerBase<MDTraits>
50
51
52
53
54
{
    using Scalar = typename MDTraits::Scalar;
    using ParentType = StaggeredCouplingManager<MDTraits>;

public:
55
56
57
58
    static constexpr auto freeFlowFaceIdx = FreeFlowPorousMediumCouplingManagerBase<MDTraits>::freeFlowFaceIdx;
    static constexpr auto freeFlowCellCenterIdx = FreeFlowPorousMediumCouplingManagerBase<MDTraits>::freeFlowCellCenterIdx;
    static constexpr auto freeFlowIdx = FreeFlowPorousMediumCouplingManagerBase<MDTraits>::freeFlowIdx;
    static constexpr auto porousMediumIdx = FreeFlowPorousMediumCouplingManagerBase<MDTraits>::porousMediumIdx;
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

private:

    using SolutionVector = typename MDTraits::SolutionVector;

    // obtain the type tags of the sub problems
    using StokesTypeTag = typename MDTraits::template SubDomain<0>::TypeTag;
    using DarcyTypeTag = typename MDTraits::template SubDomain<2>::TypeTag;

    using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >;
    using CouplingStencil = CouplingStencils::mapped_type;

    // the sub domain type tags
    template<std::size_t id>
    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;

    static constexpr bool isCompositional = GetPropType<SubDomainTypeTag<0>, Properties::ModelTraits>::numFluidComponents()> 1;

    template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
    template<std::size_t id> using NumEqVector = GetPropType<SubDomainTypeTag<id>, Properties::NumEqVector>;
    template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
    template<std::size_t id> using GridVolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>;
    template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
    template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
    template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
    template<std::size_t id> using ElementBoundaryTypes = GetPropType<SubDomainTypeTag<id>, Properties::ElementBoundaryTypes>;
    template<std::size_t id> using ElementFluxVariablesCache = typename GetPropType<SubDomainTypeTag<id>, Properties::GridFluxVariablesCache>::LocalView;
    template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
    template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
    template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
    template<std::size_t id> using PrimaryVariables = typename MDTraits::template SubDomain<id>::PrimaryVariables;
    template<std::size_t id> using SubControlVolumeFace  = typename FVElementGeometry<id>::SubControlVolumeFace;

    using CellCenterSolutionVector = GetPropType<StokesTypeTag, Properties::CellCenterSolutionVector>;

94
    using VelocityVector = typename Element<freeFlowIdx>::Geometry::GlobalCoordinate;
95
96
97
98
99
100

    using CouplingMapper = StokesDarcyCouplingMapperBox<MDTraits>;

    // Each context object contains the data related to one coupling segment
    struct StationaryStokesCouplingContext
    {
101
102
        Element<porousMediumIdx> element;
        FVElementGeometry<porousMediumIdx> fvGeometry;
103
104
        std::size_t darcyScvfIdx;
        std::size_t stokesScvfIdx;
105
106
107
        VolumeVariables<porousMediumIdx> volVars;
        std::unique_ptr< ElementVolumeVariables<porousMediumIdx> > elementVolVars;
        std::unique_ptr< ElementFluxVariablesCache<porousMediumIdx> > elementFluxVarsCache;
108
109
110
111
112
        typename CouplingMapper::CouplingSegment::Geometry segmentGeometry;
    };

    struct StationaryDarcyCouplingContext
    {
113
114
        Element<freeFlowIdx> element;
        FVElementGeometry<freeFlowIdx> fvGeometry;
115
116
117
        std::size_t stokesScvfIdx;
        std::size_t darcyScvfIdx;
        VelocityVector velocity;
118
        VolumeVariables<freeFlowIdx> volVars;
119
120
121
122
123
124
125
126
127
        typename CouplingMapper::CouplingSegment::Geometry segmentGeometry;
    };
public:

    using ParentType::couplingStencil;
    using ParentType::updateCouplingContext;
    using CouplingData = StokesDarcyCouplingData<MDTraits, StokesDarcyCouplingManager<MDTraits>>;

    //! Constructor
128
129
    StokesDarcyCouplingManagerImplementation(std::shared_ptr<const GridGeometry<freeFlowIdx>> stokesFvGridGeometry,
                                             std::shared_ptr<const GridGeometry<porousMediumIdx>> darcyFvGridGeometry)
130
131
132
133
134
135
136
137
    { }

    /*!
     * \brief Methods to be accessed by main
     */
    // \{

    //! Initialize the coupling manager
138
139
    void init(std::shared_ptr<const Problem<freeFlowIdx>> stokesProblem,
              std::shared_ptr<const Problem<porousMediumIdx>> darcyProblem,
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
              const SolutionVector& curSol)
    {
        if (Dune::FloatCmp::ne(stokesProblem->gravity(), darcyProblem->spatialParams().gravity({})))
            DUNE_THROW(Dune::InvalidStateException, "Both models must use the same gravity vector");

        this->setSubProblems(std::make_tuple(stokesProblem, stokesProblem, darcyProblem));
        this->curSol() = curSol;
        couplingData_ = std::make_shared<CouplingData>(*this);
        computeStencils();
    }

    //! Update after the grid has changed
    void update()
    { }

    // \}

    //! Update the solution vector before assembly
    void updateSolution(const SolutionVector& curSol)
    { this->curSol() = curSol; }

    //! Prepare the coupling stencils
    void computeStencils()
    {
        couplingMapper_.computeCouplingMapsAndStencils(*this,
                                                       darcyToStokesCellCenterCouplingStencils_,
                                                       darcyToStokesFaceCouplingStencils_,
                                                       stokesCellCenterCouplingStencils_,
                                                       stokesFaceCouplingStencils_);

        for(auto&& stencil : darcyToStokesCellCenterCouplingStencils_)
            removeDuplicates_(stencil.second);
        for(auto&& stencil : darcyToStokesFaceCouplingStencils_)
        {
            removeDuplicates_(stencil.second.first);
            removeDuplicates_(stencil.second.second);
        }
        for(auto&& stencil : stokesCellCenterCouplingStencils_)
            removeDuplicates_(stencil.second);
        for(auto&& stencil : stokesFaceCouplingStencils_)
            removeDuplicates_(stencil.second);
    }

    /*!
     * \brief Methods to be accessed by the assembly
     */
    // \{

    using ParentType::evalCouplingResidual;

    /*!
     * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Darcy information)
     */
193
194
    template<std::size_t i, class Assembler, std::enable_if_t<(i == freeFlowCellCenterIdx || i == freeFlowFaceIdx), int> = 0>
    void bindCouplingContext(Dune::index_constant<i> domainI, const Element<freeFlowCellCenterIdx>& element, const Assembler& assembler) const
195
196
197
    {
        stokesCouplingContext_.clear();

198
        const auto stokesElementIdx = this->problem(freeFlowIdx).gridGeometry().elementMapper().index(element);
199
200
201
202
203
204
205
206
        boundStokesElemIdx_ = stokesElementIdx;

        // do nothing if the element is not coupled to the other domain
        if(!couplingMapper_.stokesElementToDarcyElementMap().count(stokesElementIdx))
            return;

        // prepare the coupling context
        const auto& couplingSegments = couplingMapper_.stokesElementToDarcyElementMap().at(stokesElementIdx);
207
        auto darcyFvGeometry = localView(this->problem(porousMediumIdx).gridGeometry());
208
209
210

        for(const auto& couplingSegment : couplingSegments)
        {
211
            const auto& darcyElement = this->problem(porousMediumIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingSegment.eIdx);
212
213
214
215
            darcyFvGeometry.bind(darcyElement);
            const auto& scvf = darcyFvGeometry.scvf(couplingSegment.scvfIdx);
            const auto& scv = darcyFvGeometry.scv(scvf.insideScvIdx());

216
217
            auto darcyElemVolVars = localView(assembler.gridVariables(porousMediumIdx).curGridVolVars());
            auto darcyElemFluxVarsCache = localView(assembler.gridVariables(porousMediumIdx).gridFluxVarsCache());
218

219
220
            darcyElemVolVars.bind(darcyElement, darcyFvGeometry, this->curSol()[porousMediumIdx]);
            darcyElemVolVars.bindElement(darcyElement, darcyFvGeometry, this->curSol()[porousMediumIdx]);
221
222
            darcyElemFluxVarsCache.bind(darcyElement, darcyFvGeometry, darcyElemVolVars);

223
224
225
            const auto darcyElemSol = elementSolution(darcyElement, this->curSol()[porousMediumIdx], this->problem(porousMediumIdx).gridGeometry());
            VolumeVariables<porousMediumIdx> darcyVolVars;
            darcyVolVars.update(darcyElemSol, this->problem(porousMediumIdx), darcyElement, scv);
226
227
228

            // add the context
            stokesCouplingContext_.push_back({darcyElement, darcyFvGeometry, couplingSegment.scvfIdx, couplingSegment.flipScvfIdx, darcyVolVars,
229
230
                                              std::make_unique< ElementVolumeVariables<porousMediumIdx> >( std::move(darcyElemVolVars)),
                                              std::make_unique< ElementFluxVariablesCache<porousMediumIdx> >( std::move(darcyElemFluxVarsCache)),
231
232
233
234
235
236
237
238
                                              couplingSegment.geometry});
        }
    }

    /*!
     * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information)
     */
    template<class Assembler>
239
    void bindCouplingContext(Dune::index_constant<porousMediumIdx> domainI, const Element<porousMediumIdx>& element, const Assembler& assembler) const
240
241
242
    {
        darcyCouplingContext_.clear();

243
        const auto darcyElementIdx = this->problem(porousMediumIdx).gridGeometry().elementMapper().index(element);
244
245
246
247
248
249
250
251
        boundDarcyElemIdx_ = darcyElementIdx;

        // do nothing if the element is not coupled to the other domain
        if(!couplingMapper_.darcyElementToStokesElementMap().count(darcyElementIdx))
            return;

        // prepare the coupling context
        const auto& couplingSegments = couplingMapper_.darcyElementToStokesElementMap().at(darcyElementIdx);
252
        auto stokesFvGeometry = localView(this->problem(freeFlowIdx).gridGeometry());
253
254
255

        for(const auto& couplingSegment : couplingSegments)
        {
256
            const auto& stokesElement = this->problem(freeFlowIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingSegment.eIdx);
257
258
259
260
261
262
263
            stokesFvGeometry.bindElement(stokesElement);

            VelocityVector faceVelocity(0.0);

            for(const auto& scvf : scvfs(stokesFvGeometry))
            {
                if(scvf.index() == couplingSegment.scvfIdx)
264
                    faceVelocity[scvf.directionIndex()] = this->curSol()[freeFlowFaceIdx][scvf.dofIndex()];
265
266
            }

267
268
            using PriVarsType = typename VolumeVariables<freeFlowCellCenterIdx>::PrimaryVariables;
            const auto& cellCenterPriVars = this->curSol()[freeFlowCellCenterIdx][couplingSegment.eIdx];
269
270
            const auto elemSol = makeElementSolutionFromCellCenterPrivars<PriVarsType>(cellCenterPriVars);

271
            VolumeVariables<freeFlowIdx> stokesVolVars;
272
            for(const auto& scv : scvs(stokesFvGeometry))
273
                stokesVolVars.update(elemSol, this->problem(freeFlowIdx), stokesElement, scv);
274
275
276
277
278
279
280
281
282
283
284
285

            // add the context
            darcyCouplingContext_.push_back({stokesElement, stokesFvGeometry,
                                             couplingSegment.scvfIdx, couplingSegment.flipScvfIdx,
                                             faceVelocity, stokesVolVars, couplingSegment.geometry});
        }
    }

    /*!
     * \brief Update the coupling context for the Darcy residual w.r.t. Darcy DOFs
     */
    template<class LocalAssemblerI>
286
    void updateCouplingContext(Dune::index_constant<porousMediumIdx> domainI,
287
                               const LocalAssemblerI& localAssemblerI,
288
                               Dune::index_constant<porousMediumIdx> domainJ,
289
                               std::size_t dofIdxGlobalJ,
290
                               const PrimaryVariables<porousMediumIdx>& priVarsJ,
291
292
293
294
295
296
297
298
299
                               int pvIdxJ)
    {
        this->curSol()[domainJ][dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
    }

    /*!
     * \brief Update the coupling context for the Darcy residual w.r.t. the Stokes cell-center DOFs (DarcyToCC)
     */
    template<class LocalAssemblerI>
300
    void updateCouplingContext(Dune::index_constant<porousMediumIdx> domainI,
301
                               const LocalAssemblerI& localAssemblerI,
302
                               Dune::index_constant<freeFlowCellCenterIdx> domainJ,
303
                               const std::size_t dofIdxGlobalJ,
304
                               const PrimaryVariables<freeFlowCellCenterIdx>& priVars,
305
306
307
308
309
310
                               int pvIdxJ)
    {
        this->curSol()[domainJ][dofIdxGlobalJ] = priVars;

        for (auto& data : darcyCouplingContext_)
        {
311
            const auto stokesElemIdx = this->problem(freeFlowIdx).gridGeometry().elementMapper().index(data.element);
312
313
314
315

            if(stokesElemIdx != dofIdxGlobalJ)
                continue;

316
            using PriVarsType = typename VolumeVariables<freeFlowCellCenterIdx>::PrimaryVariables;
317
318
319
            const auto elemSol = makeElementSolutionFromCellCenterPrivars<PriVarsType>(priVars);

            for(const auto& scv : scvs(data.fvGeometry))
320
                data.volVars.update(elemSol, this->problem(freeFlowIdx), data.element, scv);
321
322
323
324
325
326
327
        }
    }

    /*!
     * \brief Update the coupling context for the Darcy residual w.r.t. the Stokes face DOFs (DarcyToFace)
     */
    template<class LocalAssemblerI>
328
    void updateCouplingContext(Dune::index_constant<porousMediumIdx> domainI,
329
                               const LocalAssemblerI& localAssemblerI,
330
                               Dune::index_constant<freeFlowFaceIdx> domainJ,
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
                               const std::size_t dofIdxGlobalJ,
                               const PrimaryVariables<1>& priVars,
                               int pvIdxJ)
    {
        this->curSol()[domainJ][dofIdxGlobalJ] = priVars;

        for (auto& data : darcyCouplingContext_)
        {
            for(const auto& scvf : scvfs(data.fvGeometry))
            {
                if(scvf.dofIndex() == dofIdxGlobalJ)
                    data.velocity[scvf.directionIndex()] = priVars;
            }
        }
    }

    /*!
     * \brief Update the coupling context for the Stokes cc residual w.r.t. the Darcy DOFs (FaceToDarcy)
     */
350
    template<std::size_t i, class LocalAssemblerI, std::enable_if_t<(i == freeFlowCellCenterIdx || i == freeFlowFaceIdx), int> = 0>
351
352
    void updateCouplingContext(Dune::index_constant<i> domainI,
                               const LocalAssemblerI& localAssemblerI,
353
                               Dune::index_constant<porousMediumIdx> domainJ,
354
                               const std::size_t dofIdxGlobalJ,
355
                               const PrimaryVariables<porousMediumIdx>& priVars,
356
357
358
359
360
361
362
                               int pvIdxJ)
    {
        this->curSol()[domainJ][dofIdxGlobalJ] = priVars;

        for (auto& data : stokesCouplingContext_)
        {
            //Derivatives are assumed to be always calculated with respect to unknowns associated with its own element
363
            const auto darcyElemSol = elementSolution(data.element, this->curSol()[porousMediumIdx], this->problem(porousMediumIdx).gridGeometry());
364
365
366
367

            const auto& scvf = data.fvGeometry.scvf(data.darcyScvfIdx);
            const auto& scv = data.fvGeometry.scv(scvf.insideScvIdx());
            if(scv.dofIndex() == dofIdxGlobalJ)
368
                data.volVars.update(darcyElemSol, this->problem(porousMediumIdx), data.element, scv);
369
370
371
372
373
374

            for (auto&& scv : scvs(data.fvGeometry))
            {
                if(scv.dofIndex() == dofIdxGlobalJ)
                {
                    auto& volVars = (*data.elementVolVars)[scv];
375
                    volVars.update(darcyElemSol, this->problem(porousMediumIdx), data.element, scv);
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
                }
            }
        }
    }

    // \}

    /*!
     * \brief Access the coupling data
     */
    const auto& couplingData() const
    {
        return *couplingData_;
    }

    /*!
     * \brief Access the coupling context needed for the Stokes domain
     */
394
    const auto& stokesCouplingContext(const Element<freeFlowIdx>& element, const SubControlVolumeFace<freeFlowIdx>& scvf) const
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
    {
        if (stokesCouplingContext_.empty() || boundStokesElemIdx_ != scvf.insideScvIdx())
            DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());

        for(const auto& context : stokesCouplingContext_)
        {
            if(scvf.index() == context.stokesScvfIdx)
                return context;
        }

        DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
    }

    /*!
     * \brief Access the coupling context needed for the Darcy domain
     */
411
    const auto& darcyCouplingContext(const Element<porousMediumIdx>& element, const SubControlVolumeFace<porousMediumIdx>& scvf) const
412
    {
413
        if (darcyCouplingContext_.empty() || boundDarcyElemIdx_ != this->problem(porousMediumIdx).gridGeometry().elementMapper().index(element))
414
415
416
417
418
419
420
421
422
423
424
425
426
427
            DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());

        for(const auto& context : darcyCouplingContext_)
        {
            if(scvf.index() == context.darcyScvfIdx)
                return context;
        }

        DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
    }

    /*!
     * \brief Access the coupling context needed for the Darcy domain
     */
428
    const auto& darcyCouplingContextVector(const Element<porousMediumIdx>& element, const SubControlVolumeFace<porousMediumIdx>& scvf) const
429
    {
430
        if (darcyCouplingContext_.empty() || boundDarcyElemIdx_ != this->problem(porousMediumIdx).gridGeometry().elementMapper().index(element))
431
432
433
434
435
436
437
438
            DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());

        return darcyCouplingContext_;
    }

    /*!
     * \brief Access the coupling context needed for the Stokes domain
     */
439
    const auto& stokesCouplingContextVector(const Element<freeFlowIdx>& element, const SubControlVolumeFace<freeFlowIdx>& scvf) const
440
441
442
443
444
445
446
447
448
449
450
451
452
453
    {
        if (stokesCouplingContext_.empty() || boundStokesElemIdx_ != scvf.insideScvIdx())
            DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
        return stokesCouplingContext_;
    }

    /*!
     * \brief The coupling stencils
     */
    // \{

    /*!
     * \brief The Stokes cell center coupling stencil w.r.t. Darcy DOFs
     */
454
455
456
    const CouplingStencil& couplingStencil(Dune::index_constant<freeFlowCellCenterIdx> domainI,
                                           const Element<freeFlowIdx>& element,
                                           Dune::index_constant<porousMediumIdx> domainJ) const
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
    {
        const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
        if(stokesCellCenterCouplingStencils_.count(eIdx))
            return stokesCellCenterCouplingStencils_.at(eIdx);
        else
            return emptyStencil_;
    }

    /*!
     * \brief The coupling stencil of domain I, i.e. which domain J DOFs
     *        the given domain I element's residual depends on.
     */
    template<std::size_t i, std::size_t j>
    const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI,
                                           const Element<i>& element,
                                           Dune::index_constant<j> domainJ) const
    { return emptyStencil_; }

    /*!
     * \brief The coupling stencil of domain I, i.e. which domain J dofs
     *        the given domain I element's residual depends on.
     */
479
480
481
    const CouplingStencil& couplingStencil(Dune::index_constant<porousMediumIdx> domainI,
                                           const Element<porousMediumIdx>& element,
                                           Dune::index_constant<freeFlowCellCenterIdx> domainJ) const
482
483
484
485
486
487
488
489
490
491
492
493
    {
        const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
        if(darcyToStokesCellCenterCouplingStencils_.count(eIdx))
            return darcyToStokesCellCenterCouplingStencils_.at(eIdx);
        else
            return emptyStencil_;
    }

    /*!
     * \brief The coupling stencil of domain I, i.e. which domain J dofs
     *        the given domain I element's residual depends on.
     */
494
495
496
    const CouplingStencil& couplingStencil(Dune::index_constant<porousMediumIdx> domainI,
                                           const Element<porousMediumIdx>& element,
                                           Dune::index_constant<freeFlowFaceIdx> domainJ) const
497
498
499
500
501
502
503
504
505
506
507
508
509
510
    {
        const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
        if (darcyToStokesFaceCouplingStencils_.count(eIdx))
            return darcyToStokesFaceCouplingStencils_.at(eIdx).first;
        else
            return emptyStencil_;
    }

    /*!
     * \brief The coupling stencil of domain I, i.e. which domain J DOFs
     *        the given domain I element's residual depends on.
     */
    template<std::size_t i, std::size_t j>
    const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI,
511
                                           const SubControlVolumeFace<freeFlowIdx>& scvf,
512
513
514
515
516
517
                                           Dune::index_constant<j> domainJ) const
    { return emptyStencil_; }

    /*!
     * \brief The coupling stencil of a Stokes face w.r.t. Darcy DOFs
     */
518
519
520
    const CouplingStencil& couplingStencil(Dune::index_constant<freeFlowFaceIdx> domainI,
                                           const SubControlVolumeFace<freeFlowIdx>& scvf,
                                           Dune::index_constant<porousMediumIdx> domainJ) const
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
    {
        const auto faceDofIdx = scvf.dofIndex();
        if(stokesFaceCouplingStencils_.count(faceDofIdx))
            return stokesFaceCouplingStencils_.at(faceDofIdx);
        else
            return emptyStencil_;
    }

    // \}

    /*!
     * \brief There are no additional dof dependencies
     */
    template<class IdType>
    const std::vector<std::size_t>& getAdditionalDofDependencies(IdType id, std::size_t stokesElementIdx) const
    { return emptyStencil_; }

    /*!
     * \brief There are no additional dof dependencies
     */
    template<class IdType>
    const std::vector<std::size_t>& getAdditionalDofDependenciesInverse(IdType id, std::size_t darcyElementIdx) const
    { return emptyStencil_; }

    /*!
     * \brief Returns whether a given free flow scvf is coupled to the other domain
     */
548
    bool isCoupledEntity(Dune::index_constant<freeFlowIdx>, const SubControlVolumeFace<freeFlowFaceIdx>& scvf) const
549
550
551
552
553
554
555
    {
        return stokesFaceCouplingStencils_.count(scvf.dofIndex());
    }

    /*!
     * \brief Returns whether a given free flow scvf is coupled to the other domain
     */
556
    bool isCoupledEntity(Dune::index_constant<porousMediumIdx>, const Element<porousMediumIdx>& element, const SubControlVolumeFace<porousMediumIdx>& scvf) const
557
    {
558
        return couplingMapper_.isCoupledDarcyScvf(this->problem(porousMediumIdx).gridGeometry().elementMapper().index(element),
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
                                                  scvf.index());
    }

    //! compute the shape function for a given point and geometry
    template<class FVGG, class Geometry, class GlobalPosition, class ShapeValues>
    void getShapeValues(const FVGG& gridGeometry, const Geometry& geo, const GlobalPosition& globalPos, ShapeValues& shapeValues) const
    {
        const auto ipLocal = geo.local(globalPos);
        const auto& localBasis = gridGeometry.feCache().get(geo.type()).localBasis();
        localBasis.evaluateFunction(ipLocal, shapeValues);
    }

protected:

    //! Return a reference to an empty stencil
    std::vector<std::size_t>& emptyStencil()
    { return emptyStencil_; }

    void removeDuplicates_(std::vector<std::size_t>& stencil)
    {
        std::sort(stencil.begin(), stencil.end());
        stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end());
    }

private:

    std::vector<bool> isCoupledDarcyDof_;
    std::shared_ptr<CouplingData> couplingData_;

    std::unordered_map<std::size_t, std::vector<std::size_t> > stokesCellCenterCouplingStencils_;
    std::unordered_map<std::size_t, std::vector<std::size_t> > stokesFaceCouplingStencils_;
    std::unordered_map<std::size_t, std::vector<std::size_t> > darcyToStokesCellCenterCouplingStencils_;
    std::unordered_map<std::size_t, std::pair<std::vector<std::size_t>,std::vector<std::size_t>> > darcyToStokesFaceCouplingStencils_;
    std::vector<std::size_t> emptyStencil_;

    ////////////////////////////////////////////////////////////////////////////
    //! The coupling context
    ////////////////////////////////////////////////////////////////////////////
    mutable std::vector<StationaryStokesCouplingContext> stokesCouplingContext_;
    mutable std::vector<StationaryDarcyCouplingContext> darcyCouplingContext_;

    mutable std::size_t boundStokesElemIdx_;
    mutable std::size_t boundDarcyElemIdx_;

    CouplingMapper couplingMapper_;
};

} // end namespace Dumux

#endif