fvelementgeometry.hh 15.5 KB
Newer Older
1
2
3
4
5
6
7
// -*- 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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
13
14
15
16
17
18
19
20
 *   (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
21
 * \ingroup FacetCoupling
22
23
24
25
26
27
 * \brief \copydoc Dumux::BoxFacetCouplingFVElementGeometry
 */
#ifndef DUMUX_FACETCOUPLING_BOX_FV_ELEMENT_GEOMETRY_HH
#define DUMUX_FACETCOUPLING_BOX_FV_ELEMENT_GEOMETRY_HH

#include <algorithm>
28
#include <optional>
29

30
#include <dune/geometry/type.hh>
31

32
#include <dumux/common/indextraits.hh>
33
34
35
36
37
38
#include <dumux/discretization/scvandscvfiterators.hh>
#include <dumux/discretization/box/boxgeometryhelper.hh>

namespace Dumux {

/*!
39
 * \ingroup FacetCoupling
40
41
42
43
44
 * \brief Base class for the element-local finite volume geometry for box models
 *        in the context of models considering coupling of different domains across the
 *        bulk grid facets. This builds up the sub control volumes and sub control volume
 *        faces for an element.
 * \tparam GG the finite volume grid geometry type
45
 * \tparam enableGridGeometryCache if the grid geometry is cached or not
46
 */
47
template<class GG, bool enableGridGeometryCache>
48
49
50
51
52
53
54
55
56
class BoxFacetCouplingFVElementGeometry;

//! specialization in case the FVElementGeometries are stored
template<class GG>
class BoxFacetCouplingFVElementGeometry<GG, true>
{
    using GridView = typename GG::GridView;
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
57
58
    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
    using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
59
60
61
    using CoordScalar = typename GridView::ctype;
    using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
public:
62
63
    //! export type of the element
    using Element = typename GridView::template Codim<0>::Entity;
64
65
66
67
68
    //! export type of subcontrol volume
    using SubControlVolume = typename GG::SubControlVolume;
    //! export type of subcontrol volume face
    using SubControlVolumeFace = typename GG::SubControlVolumeFace;
    //! export type of finite volume grid geometry
69
    using GridGeometry = GG;
70
71
72
73
    //! the maximum number of scvs per element (2^dim for cubes)
    static constexpr std::size_t maxNumElementScvs = (1<<dim);

    //! Constructor
74
75
    BoxFacetCouplingFVElementGeometry(const GridGeometry& gridGeometry)
    : gridGeometryPtr_(&gridGeometry)
76
    {}
77
78

    //! Get a sub control volume with a local scv index
79
    const SubControlVolume& scv(LocalIndexType scvIdx) const
80
    { return gridGeometry().scvs(eIdx_)[scvIdx]; }
81
82

    //! Get a sub control volume face with a local scvf index
83
    const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const
84
    { return gridGeometry().scvfs(eIdx_)[scvfIdx]; }
85
86
87
88
89
90
91
92
93

    //! iterator range for sub control volumes. Iterates over
    //! all scvs of the bound element.
    //! This is a free function found by means of ADL
    //! To iterate over all sub control volumes of this FVElementGeometry use
    //! for (auto&& scv : scvs(fvGeometry))
    friend inline Dune::IteratorRange<typename std::vector<SubControlVolume>::const_iterator>
    scvs(const BoxFacetCouplingFVElementGeometry& fvGeometry)
    {
94
        const auto& g = fvGeometry.gridGeometry();
95
96
97
98
99
100
101
102
103
104
105
106
        using Iter = typename std::vector<SubControlVolume>::const_iterator;
        return Dune::IteratorRange<Iter>(g.scvs(fvGeometry.eIdx_).begin(), g.scvs(fvGeometry.eIdx_).end());
    }

    //! iterator range for sub control volumes faces. Iterates over
    //! all scvfs of the bound element.
    //! This is a free function found by means of ADL
    //! To iterate over all sub control volume faces of this FVElementGeometry use
    //! for (auto&& scvf : scvfs(fvGeometry))
    friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
    scvfs(const BoxFacetCouplingFVElementGeometry& fvGeometry)
    {
107
        const auto& g = fvGeometry.gridGeometry();
108
109
110
111
112
113
        using Iter = typename std::vector<SubControlVolumeFace>::const_iterator;
        return Dune::IteratorRange<Iter>(g.scvfs(fvGeometry.eIdx_).begin(), g.scvfs(fvGeometry.eIdx_).end());
    }

    //! Get a local finite element basis
    const FeLocalBasis& feLocalBasis() const
114
    { return gridGeometry().feCache().get(element_->type()).localBasis(); }
115
116
117

    //! The total number of sub control volumes
    std::size_t numScv() const
118
    { return gridGeometry().scvs(eIdx_).size(); }
119
120
121

    //! The total number of sub control volume faces
    std::size_t numScvf() const
122
    { return gridGeometry().scvfs(eIdx_).size(); }
123
124
125
126
127
128
129
130
131
132
133
134
135
136

    //! this function is for compatibility reasons with cc methods
    //! The box stencil is always element-local so bind and bindElement
    //! are identical.
    void bind(const Element& element)
    {
        this->bindElement(element);
    }

    //! Binding of an element, has to be called before using the fvgeometries
    //! Prepares all the volume variables within the element
    //! For compatibility reasons with the FVGeometry cache being disabled
    void bindElement(const Element& element)
    {
137
        element_ = element;
138
        eIdx_ = gridGeometry().elementMapper().index(element);
139
140
    }

141
142
    //! The global finite volume geometry we are a restriction of
    const GridGeometry& gridGeometry() const
143
    { return *gridGeometryPtr_; }
144

145
146
147
148
149
150
151
152
    //! Returns true if bind/bindElement has already been called
    bool isBound() const
    { return static_cast<bool>(element_); }

    //! The bound element
    const Element& element() const
    { return *element_; }

153
private:
154
    const GridGeometry* gridGeometryPtr_;
155

156
    GridIndexType eIdx_;
157
    std::optional<Element> element_;
158
159
160
161
162
163
164
165
166
167
};

//! specialization in case the geometries are not stored grid-wide
template<class GG>
class BoxFacetCouplingFVElementGeometry<GG, false>
{
    using GridView = typename GG::GridView;
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;

168
169
    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
    using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
170
171
172
173
174
175
176
177

    using CoordScalar = typename GridView::ctype;
    using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;

    using GeometryHelper = BoxGeometryHelper<GridView, dim,
                                             typename GG::SubControlVolume,
                                             typename GG::SubControlVolumeFace>;
public:
178
179
    //! export type of the element
    using Element = typename GridView::template Codim<0>::Entity;
180
181
182
183
184
    //! export type of subcontrol volume
    using SubControlVolume = typename GG::SubControlVolume;
    //! export type of subcontrol volume face
    using SubControlVolumeFace = typename GG::SubControlVolumeFace;
    //! export type of finite volume grid geometry
185
    using GridGeometry = GG;
186
187
188
189
    //! the maximum number of scvs per element (2^dim for cubes)
    static constexpr std::size_t maxNumElementScvs = (1<<dim);

    //! Constructor
190
191
    BoxFacetCouplingFVElementGeometry(const GridGeometry& gridGeometry)
    : gridGeometryPtr_(&gridGeometry)
192
    {}
193
194

    //! Get a sub control volume with a local scv index
195
    const SubControlVolume& scv(LocalIndexType scvIdx) const
196
197
198
    { return scvs_[scvIdx]; }

    //! Get a sub control volume face with a local scvf index
199
    const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
    { return scvfs_[scvfIdx]; }

    //! iterator range for sub control volumes. Iterates over
    //! all scvs of the bound element.
    //! This is a free function found by means of ADL
    //! To iterate over all sub control volumes of this FVElementGeometry use
    //! for (auto&& scv : scvs(fvGeometry))
    friend inline Dune::IteratorRange<typename std::vector<SubControlVolume>::const_iterator>
    scvs(const BoxFacetCouplingFVElementGeometry& fvGeometry)
    {
        using Iter = typename std::vector<SubControlVolume>::const_iterator;
        return Dune::IteratorRange<Iter>(fvGeometry.scvs_.begin(), fvGeometry.scvs_.end());
    }

    //! iterator range for sub control volumes faces. Iterates over
    //! all scvfs of the bound element.
    //! This is a free function found by means of ADL
    //! To iterate over all sub control volume faces of this FVElementGeometry use
    //! for (auto&& scvf : scvfs(fvGeometry))
    friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
    scvfs(const BoxFacetCouplingFVElementGeometry& fvGeometry)
    {
        using Iter = typename std::vector<SubControlVolumeFace>::const_iterator;
        return Dune::IteratorRange<Iter>(fvGeometry.scvfs_.begin(), fvGeometry.scvfs_.end());
    }

    //! Get a local finite element basis
    const FeLocalBasis& feLocalBasis() const
228
    { return gridGeometry().feCache().get(element_->type()).localBasis(); }
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250

    //! The total number of sub control volumes
    std::size_t numScv() const
    { return scvs_.size(); }

    //! The total number of sub control volume faces
    std::size_t numScvf() const
    { return scvfs_.size(); }

    //! this function is for compatibility reasons with cc methods
    //! The box stencil is always element-local so bind and bindElement
    //! are identical.
    void bind(const Element& element)
    {
        this->bindElement(element);
    }

    //! Binding of an element, has to be called before using the fvgeometries
    //! Prepares all the volume variables within the element
    //! For compatibility reasons with the FVGeometry cache being disabled
    void bindElement(const Element& element)
    {
251
        element_ = element;
252
        eIdx_ = gridGeometry().elementMapper().index(element);
253
        makeElementGeometries_();
254
255
    }

256
257
    //! The global finite volume geometry we are a restriction of
    const GridGeometry& gridGeometry() const
258
    { return *gridGeometryPtr_; }
259

260
261
262
263
264
265
266
267
    //! Returns true if bind/bindElement has already been called
    bool isBound() const
    { return static_cast<bool>(element_); }

    //! The bound element
    const Element& element() const
    { return *element_; }

268
269
private:

270
    void makeElementGeometries_()
271
272
    {
        // get the element geometry
273
274
        const auto& element = *element_;
        const auto elementGeometry = element.geometry();
275
        const auto refElement = referenceElement(elementGeometry);
276
277
278
279
280
281
282
283
284
285
286

        // get the sub control volume geometries of this element
        GeometryHelper geometryHelper(elementGeometry);

        // construct the sub control volumes
        scvs_.clear();
        scvs_.reserve(elementGeometry.corners());
        using LocalIndexType = typename SubControlVolumeFace::Traits::LocalIndexType;
        for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
            scvs_.emplace_back(geometryHelper,
                               scvLocalIdx,
287
                               eIdx_,
288
                               gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim));
289
290
291
292
293
294
295
296
297
298

        // construct the sub control volume faces
        const auto numInnerScvf = (dim==1) ? 1 : element.subEntities(dim-1);
        scvfs_.clear();
        scvfs_.reserve(numInnerScvf);

        unsigned int scvfLocalIdx = 0;
        for (; scvfLocalIdx < numInnerScvf; ++scvfLocalIdx)
        {
            // find the local scv indices this scvf is connected to
299
300
            std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
                                                         static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
301
302
303
304
305
306
307
308
309
310

            // create the sub-control volume face
            scvfs_.emplace_back(geometryHelper,
                                element,
                                elementGeometry,
                                scvfLocalIdx,
                                std::move(localScvIndices));
        }

        // construct the sub control volume faces on the domain/interior boundaries
311
312
        // skip handled facets (necessary for e.g. Dune::FoamGrid)
        std::vector<unsigned int> handledFacets;
313
        for (const auto& intersection : intersections(gridGeometry().gridView(), element))
314
        {
315
316
317
318
319
            if (std::count(handledFacets.begin(), handledFacets.end(), intersection.indexInInside()))
                continue;

            handledFacets.push_back(intersection.indexInInside());

320
321
322
323
324
325
326
327
            // determine if all corners live on the facet grid
            const auto isGeometry = intersection.geometry();
            const auto numFaceCorners = isGeometry.corners();
            const auto idxInInside = intersection.indexInInside();
            const auto boundary = intersection.boundary();

            std::vector<LocalIndexType> vIndicesLocal(numFaceCorners);
            for (int i = 0; i < numFaceCorners; ++i)
328
                vIndicesLocal[i] = static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, i, dim));
329
330

            // if all vertices are living on the facet grid, this is an interiour boundary
331
            const bool isOnFacet = gridGeometry().isOnInteriorBoundary(element, intersection);
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357

            if (isOnFacet || boundary)
            {
                for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < isGeometry.corners(); ++isScvfLocalIdx)
                {
                    // find the inside scv this scvf is belonging to (localIdx = element local vertex index)
                    std::vector<LocalIndexType> localScvIndices = {vIndicesLocal[isScvfLocalIdx], vIndicesLocal[isScvfLocalIdx]};

                    // create the sub-control volume face
                    scvfs_.emplace_back(geometryHelper,
                                        intersection,
                                        isGeometry,
                                        isScvfLocalIdx,
                                        scvfLocalIdx,
                                        std::move(localScvIndices),
                                        boundary,
                                        isOnFacet);

                    // increment local counter
                    scvfLocalIdx++;
                }
            }
        }
    }

    //! The global geometry this is a restriction of
358
    const GridGeometry* gridGeometryPtr_;
359

360
361
362
363
    //! The bound element
    GridIndexType eIdx_;
    std::optional<Element> element_;

364
365
366
367
368
369
370
371
    //! vectors to store the geometries locally after binding an element
    std::vector<SubControlVolume> scvs_;
    std::vector<SubControlVolumeFace> scvfs_;
};

} // end namespace Dumux

#endif