couplingmapper.hh 17.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
// -*- 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::StokesDarcyCouplingMapper
 */

25
26
27
28
#ifndef DUMUX_STOKES_DARCY_COUPLINGMAPPER_BOX_HH
#define DUMUX_STOKES_DARCY_COUPLINGMAPPER_BOX_HH

#include <type_traits>
29
30
31
32
33
34
35
36
#include <unordered_map>
#include <vector>

#include <dune/common/exceptions.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/geometry/affinegeometry.hh>

#include <dumux/common/properties.hh>
37
38
#include <dumux/geometry/geometryintersection.hh>
#include <dumux/geometry/intersectingentities.hh>
39
40

#include <dumux/discretization/method.hh>
41
#include <dumux/freeflow/slipcondition.hh>
42
43
44

namespace Dumux {

45
template<class CouplingFacet, class Container, class IndexSet, class Iterator>
46
47
class FacetIterator : public Dune::ForwardIteratorFacade<FacetIterator<CouplingFacet,
                                                                       Container,
48
49
                                                                       IndexSet,
                                                                       Iterator>,
50
51
                                                         const CouplingFacet>
{
52
    using ThisType = FacetIterator<CouplingFacet, Container, IndexSet, Iterator>;
53
public:
54
55
    FacetIterator(const Iterator& it, const Container& container, const IndexSet& indices, std::size_t localIdx)
    : it_(it), container_(&container), indices_(&indices), localIdx_(localIdx) {}
56

57
    FacetIterator() : it_(Iterator()), container_(nullptr), indices_(nullptr) {}
58
59
60
61
62
63
64
65
66
67
68
69

    //! dereferencing yields a coupling facet
    const CouplingFacet& dereference() const
    {
        return container_->at(*it_);
    }

    bool equals(const ThisType& other) const
    {
        return it_ == other.it_;
    }

70
71
    void increment()
    {
72
73
74
75
76
        // Case where the index container contains the indices
        if constexpr (std::is_integral_v<typename IndexSet::value_type>)
            it_++;
        // Case where the index container contains sub-containers which store the indices
        else
77
        {
78
79
80
81
82
83
84
85
86
87
88
            it_++;
            if(it_ == (*indices_)[localIdx_].end())
            {
                auto it = std::find_if(indices_->begin()+localIdx_+1, indices_->end(), [] (const auto& idx) { return idx.size() > 0; });

                if(it != indices_->end())
                {
                    it_ = it->begin();
                    localIdx_ = std::distance(indices_->begin(), it);
                }
            }
89
90
91
        }
    }

92
93
94
private:
    Iterator it_;
    const Container* container_;
95
96
    const IndexSet* indices_;
    std::size_t localIdx_{0};
97
98
};

99
100
101
102
103
104
105
106
107
108
109
110
111
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; });

112
113
    return Dune::IteratorRange<FacetIterator>(FacetIterator(itBegin->begin(), mapper.couplingFacets(), indexSet, std::distance(indexSet.begin(), itBegin)),
                                              FacetIterator(itEnd->end(), mapper.couplingFacets(), indexSet, std::distance(itEnd, std::reverse_iterator(indexSet.begin()))));
114
115
116
117
118
119
120
121
}

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>;
122
123
124
125
    const auto& indexSet = mapper.couplingFacetIdxMap(domainI).at(eIdx);
    const auto& indexSetScvf = indexSet[localScvfIdx];
    return Dune::IteratorRange<FacetIterator>(FacetIterator(indexSetScvf.begin(), mapper.couplingFacets(), indexSetScvf, localScvfIdx),
                                              FacetIterator(indexSetScvf.end(), mapper.couplingFacets(), indexSetScvf, localScvfIdx));
126
}
127

128
129
/*!
 * \ingroup StokesDarcyCoupling
130
 * \brief Coupling mapper for Stokes and Darcy domains with equal dimension when using the Box scheme for the Darcy domain.
131
132
133
134
135
 */
template<class MDTraits>
class StokesDarcyCouplingMapperBox
{
private:
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;

    template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
    template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
    template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
    template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
    template<std::size_t id> using Scvf = typename GridGeometry<id>::SubControlVolumeFace::Traits::Geometry;

    static constexpr auto freeFlowIdx = FreeFlowPorousMediumCouplingManagerBase<MDTraits>::freeFlowIdx;
    static constexpr auto porousMediumIdx = FreeFlowPorousMediumCouplingManagerBase<MDTraits>::porousMediumIdx;

    static constexpr int dim = GridView<porousMediumIdx>::dimension;
    static constexpr int dimWorld = GridView<porousMediumIdx>::dimensionworld;

    using ctype = typename Dune::PromotionTraits<typename GridView<freeFlowIdx>::ctype,
                                                 typename GridView<porousMediumIdx>::ctype>::PromotedType;

    static_assert(GridView<freeFlowIdx>::dimension == dim, "The grids must have the same dimension");
    static_assert(GridView<freeFlowIdx>::dimensionworld == dimWorld, "The grids must have the same world dimension");
    static_assert(GridGeometry<freeFlowIdx>::discMethod == DiscretizationMethod::staggered, "The free flow domain must use the staggered discretization");
    static_assert(GridGeometry<porousMediumIdx>::discMethod == DiscretizationMethod::box, "The Darcy domain must use the Box discretization");
157
158
159

public:
    // export the type describing a coupling segment
160
    struct CouplingFacet
161
162
163
164
    {
        // each intersection segment is described by a simplex geometry of codimension one
        using Geometry = Dune::AffineGeometry<ctype, dim-1, dimWorld>;

165
166
167
168
169
        std::size_t ffEIdx;
        std::size_t pmEIdx;
        std::size_t ffScvfIdx;
        std::size_t pmScvfIdx;
        std::size_t idx;
170
171
172
        Geometry geometry;
    };

173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
    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));
    // }
197

198
199
200
    /*!
     * \brief Main update routine
     */
201
    template<class CouplingManager, class StencilA, class StencilB>
202
    void computeCouplingMapsAndStencils(const CouplingManager& couplingManager,
203
204
205
206
                                        StencilA& darcyToStokesCellCenterStencils,
                                        StencilB& darcyToStokesFaceStencils,
                                        StencilA& stokesCellCenterToDarcyStencils,
                                        StencilA& stokesFaceToDarcyStencils)
207
208
209
    {
        computeCouplingMaps(couplingManager);

210
211
        const auto& stokesProblem = couplingManager.problem(CouplingManager::freeFlowIdx);
        const auto& darcyProblem = couplingManager.problem(CouplingManager::porousMediumIdx);
212
213
214
215

        const auto& stokesFvGridGeometry = stokesProblem.gridGeometry();
        const auto& darcyFvGridGeometry = darcyProblem.gridGeometry();

216
        for(const auto& elementData : ffCouplingFacetIdxMap_)
217
        {
218
            const auto stokesEIdx = elementData.first;
219

220
            for(const auto& couplingFacet : Dumux::couplingFacets(Dune::index_constant<freeFlowIdx>(), *this, stokesEIdx))
221
            {
222
223
224
                const auto darcyEIdx = couplingFacet.pmEIdx;
                const auto stokesScvfIdx = couplingFacet.ffScvfIdx;
                const auto& stokesScvf = stokesFvGridGeometry.scvf(stokesScvfIdx);
225

226
227
228
                const auto& darcyElement = darcyFvGridGeometry.element(darcyEIdx);
                auto darcyFvGeometry = localView(darcyFvGridGeometry);
                darcyFvGeometry.bind(darcyElement);
229

230
231
232
                const auto stokesElement = stokesFvGridGeometry.element(stokesEIdx);
                auto stokesFvGeometry = localView(stokesFvGridGeometry);
                stokesFvGeometry.bind(stokesElement);
233

234
235
236
                darcyToStokesCellCenterStencils[darcyEIdx].push_back(stokesEIdx);
                darcyToStokesFaceStencils[darcyEIdx].first.push_back(stokesScvf.dofIndex());
                darcyToStokesFaceStencils[darcyEIdx].second.push_back(stokesScvf.index());
237

238
239
240
241
242
                for (auto&& scv : scvs(darcyFvGeometry))
                {
                    stokesCellCenterToDarcyStencils[stokesEIdx].push_back(scv.dofIndex());
                    stokesFaceToDarcyStencils[stokesScvf.dofIndex()].push_back(scv.dofIndex());
                }
243

244
245
246
247
248
                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)
249
                    {
250
251
252
253
254
                        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());
255
                    }
256
                }
257
258
259
260
261
262
263
            }
        }
    }

    template<class CouplingManager>
    void computeCouplingMaps(const CouplingManager& couplingManager)
    {
264
265
        const auto& stokesProblem = couplingManager.problem(CouplingManager::freeFlowIdx);
        const auto& darcyProblem = couplingManager.problem(CouplingManager::porousMediumIdx);
266
267
268
269

        const auto& stokesFvGridGeometry = stokesProblem.gridGeometry();
        const auto& darcyFvGridGeometry = darcyProblem.gridGeometry();

270
        std::size_t couplingFaceIdx = 0;
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
        // find all darcy faces coupling to stokes
        for (const auto& darcyElement : elements(darcyFvGridGeometry.gridView()))
        {
            const auto darcyEIdx = darcyFvGridGeometry.elementMapper().index(darcyElement);
            auto darcyFvGeometry = localView(darcyFvGridGeometry);
            darcyFvGeometry.bindElement(darcyElement);

            for (const auto& darcyScvf : scvfs(darcyFvGeometry))
            {
                if (!darcyScvf.boundary())
                    continue;

                // find all stokes elements that intersect with the face
                const auto& darcyScvfGeometry = darcyScvf.geometry();
                const auto rawIntersections = intersectingEntities(darcyScvfGeometry, stokesFvGridGeometry.boundingBoxTree());
                if (rawIntersections.empty())
                    continue;

                for (const auto& rawIntersection : rawIntersections)
                {
                    const auto stokesEIdx = rawIntersection.second();
                    const auto stokesElement = stokesFvGridGeometry.element(stokesEIdx);
                    auto stokesFvGeometry = localView(stokesFvGridGeometry);
                    stokesFvGeometry.bindElement(stokesElement);

                    for (const auto& stokesScvf : scvfs(stokesFvGeometry))
                    {
                        if (!stokesScvf.boundary())
                            continue;

                        // intersect the geometries
302
                        using IntersectionAlgorithm = GeometryIntersection<Scvf<porousMediumIdx>, Scvf<freeFlowIdx>>;
303
304
305
                        typename IntersectionAlgorithm::Intersection rawIs;
                        if(IntersectionAlgorithm::intersection(darcyScvfGeometry, stokesScvf.geometry(), rawIs))
                        {
306
307
308
309
310
311
312
313
314
315
316
                            if(pmCouplingFacetIdxMap_[darcyEIdx].size() == 0)
                                pmCouplingFacetIdxMap_[darcyEIdx].resize(darcyFvGeometry.numScvf());

                            if(ffCouplingFacetIdxMap_[stokesEIdx].size() == 0)
                                ffCouplingFacetIdxMap_[stokesEIdx].resize(stokesFvGeometry.numScvf());

                            const auto is = typename CouplingFacet::Geometry(Dune::GeometryTypes::simplex(dim-1), rawIs);
                            pmCouplingFacetIdxMap_[darcyEIdx][darcyScvf.index()].push_back(couplingFaceIdx);
                            ffCouplingFacetIdxMap_[stokesEIdx][stokesScvf.localFaceIdx()].push_back(couplingFaceIdx);
                            couplingFacets_.push_back({stokesEIdx, darcyEIdx, stokesScvf.index(), darcyScvf.index(), couplingFaceIdx, is});
                            couplingFaceIdx++;
317
318
319
320
321
322
323
324
325
326
327
328
                        }
                    }
                }
            }
        }
    }

    /*!
     * \brief Returns whether a Darcy scvf is coupled to the other domain
     */
    bool isCoupledDarcyScvf(std::size_t eIdx, std::size_t scvfLocalIdx) const
    {
329
330
        return (pmCouplingFacetIdxMap_.count(eIdx) == 0) ?  false
                                                         : pmCouplingFacetIdxMap_.at(eIdx)[scvfLocalIdx].size() > 0;
331
332
333
    }

    /*!
334
     * \brief A map that returns all elements coupled to the other domain
335
     */
336
337
    template<std::size_t i>
    const auto& couplingFacetIdxMap(Dune::index_constant<i> domainI) const
338
    {
339
340
341
342
        if constexpr (i == porousMediumIdx)
            return pmCouplingFacetIdxMap_;
        else if constexpr (i == freeFlowIdx)
            return ffCouplingFacetIdxMap_;
343
344
    }

345
    const CouplingFacet& couplingFacet(std::size_t idx) const
346
    {
347
        return couplingFacets_[idx];
348
349
    }

350
    const FacetContainer& couplingFacets() const
351
    {
352
        return couplingFacets_;
353
    }
354

355

356
private:
357
358
359
    FacetContainer couplingFacets_;
    std::unordered_map<std::size_t, IndexContainerElement> pmCouplingFacetIdxMap_;
    std::unordered_map<std::size_t, IndexContainerElement> ffCouplingFacetIdxMap_;
360
361
362
363
364
};

} // end namespace Dumux

#endif