interactionvolume.hh 15.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- 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 2 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
Dennis Gläser's avatar
Dennis Gläser committed
21
 * \ingroup CCMpfaDiscretization
22
 * \brief Class for the interaction volume of the mpfa-o scheme.
23
24
25
26
 */
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH

Dennis Gläser's avatar
Dennis Gläser committed
27
28
29
30
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
#include <dune/common/fvector.hh>

31
#include <dumux/common/math.hh>
Dennis Gläser's avatar
Dennis Gläser committed
32
#include <dumux/common/properties.hh>
33

Dennis Gläser's avatar
Dennis Gläser committed
34
#include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh>
35
#include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh>
Dennis Gläser's avatar
Dennis Gläser committed
36
#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
Dennis Gläser's avatar
Dennis Gläser committed
37
38
#include <dumux/discretization/cellcentered/mpfa/localfacedata.hh>
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
39
40

#include "localsubcontrolentities.hh"
Dennis Gläser's avatar
Dennis Gläser committed
41
#include "interactionvolumeindexset.hh"
42
43
44

namespace Dumux
{
Dennis Gläser's avatar
Dennis Gläser committed
45
46
47
48
//! Forward declaration of the interaction volume specialization for the mpfa-o scheme
template< class TypeTag >
class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>;

49
//! Specialization of the interaction volume traits class for the mpfa-o method
Dennis Gläser's avatar
Dennis Gläser committed
50
51
template< class TypeTag >
struct CCMpfaInteractionVolumeTraits< CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> >
52
{
Dennis Gläser's avatar
Dennis Gläser committed
53
54
55
private:
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
56

Dennis Gläser's avatar
Dennis Gläser committed
57
58
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
59

Dennis Gläser's avatar
Dennis Gläser committed
60
    using InteractionVolumeType = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>;
Dennis Gläser's avatar
Dennis Gläser committed
61

Dennis Gläser's avatar
Dennis Gläser committed
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
public:
    //! Per default, we use the same interaction volume everywhere (also on boundaries etc...)
    using SecondaryInteractionVolume = InteractionVolumeType;

    //! export the type used for local indices
    using LocalIndexType = std::uint8_t;
    //! export the type used for indices on the grid
    using GridIndexType = typename GridView::IndexSet::IndexType;
    //! export the type for the interaction volume index set
    using IndexSet = CCMpfaOInteractionVolumeIndexSet< DualGridNodalIndexSet<GridIndexType, LocalIndexType, dim> >;
    //! export the type used for global coordinates
    using GlobalPosition = Dune::FieldVector< typename GridView::ctype, dimWorld >;
    //! export the type of interaction-volume local scvs
    using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, GlobalPosition, dim >;
    //! export the type of interaction-volume local scvfs
    using LocalScvfType = CCMpfaOInteractionVolumeLocalScvf< IndexSet >;
    //! export the type of used for the iv-local face data
    using LocalFaceData = InteractionVolumeLocalFaceData<GridIndexType, LocalIndexType>;
    //! export the type of face data container (use dynamic container here)
    using LocalFaceDataContainer = std::vector< LocalFaceData >;
    //! export the type used for iv-local matrices
    using Matrix = Dune::DynamicMatrix< Scalar >;
    //! export the type used for iv-local vectors
    using Vector = Dune::DynamicVector< Scalar >;
86
87
    //! export the type used for the iv-stencils
    using Stencil = std::vector< GridIndexType >;
Dennis Gläser's avatar
Dennis Gläser committed
88
89
    //! export the data handle type for this iv
    using DataHandle = InteractionVolumeDataHandle< TypeTag, InteractionVolumeType >;
90
};
91

92
/*!
Dennis Gläser's avatar
Dennis Gläser committed
93
94
 * \ingroup CCMpfaDiscretization
 * \brief Class for the interaction volume of the mpfa-o method.
95
96
 */
template<class TypeTag>
Dennis Gläser's avatar
Dennis Gläser committed
97
98
99
class CCMpfaInteractionVolumeImplementation< TypeTag, MpfaMethods::oMethod >
           : public CCMpfaInteractionVolumeBase< TypeTag,
                                                 CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> >
100
{
Dennis Gläser's avatar
Dennis Gläser committed
101
102
    using ThisType = CCMpfaInteractionVolumeImplementation< TypeTag, MpfaMethods::oMethod >;
    using ParentType = CCMpfaInteractionVolumeBase< TypeTag, ThisType >;
103

104
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
Dennis Gläser's avatar
Dennis Gläser committed
105
    using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
106
107
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Element = typename GridView::template Codim<0>::Entity;
Dennis Gläser's avatar
Dennis Gläser committed
108
109

    static constexpr int dim = GridView::dimension;
110
    using DimVector = Dune::FieldVector<Scalar, dim>;
111

Dennis Gläser's avatar
Dennis Gläser committed
112
113
114
115
116
    using TraitsType = CCMpfaInteractionVolumeTraits< ThisType >;
    using Matrix = typename TraitsType::Matrix;
    using LocalScvType = typename TraitsType::LocalScvType;
    using LocalScvfType = typename TraitsType::LocalScvfType;
    using LocalFaceData = typename TraitsType::LocalFaceData;
117

Dennis Gläser's avatar
Dennis Gläser committed
118
119
120
    using IndexSet = typename TraitsType::IndexSet;
    using GridIndexType = typename TraitsType::GridIndexType;
    using LocalIndexType = typename TraitsType::LocalIndexType;
121
    using Stencil = typename TraitsType::Stencil;
122

Dennis Gläser's avatar
Dennis Gläser committed
123
124
125
126
127
128
129
130
131
132
133
134
135
    //! Data attached to scvf touching Dirichlet boundaries.
    //! For the default o-scheme, we only store the corresponding vol vars index.
    class DirichletData
    {
        GridIndexType volVarIndex_;

    public:
        //! Constructor
        DirichletData(const GridIndexType index) : volVarIndex_(index) {}

        //! Return corresponding vol var index
        GridIndexType volVarIndex() const { return volVarIndex_; }
    };
Dennis Gläser's avatar
Dennis Gläser committed
136

Dennis Gläser's avatar
Dennis Gläser committed
137
public:
138
139
140
    //! publicly state the mpfa-scheme this interaction volume is associated with
    static constexpr MpfaMethods MpfaMethod = MpfaMethods::oMethod;

Dennis Gläser's avatar
Dennis Gläser committed
141
142
143
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
Dennis Gläser's avatar
Dennis Gläser committed
144

Dennis Gläser's avatar
Dennis Gläser committed
145
    //! Sets up the local scope for a given iv index set
146
147
148
    void setUpLocalScope(const IndexSet& indexSet,
                         const Problem& problem,
                         const FVElementGeometry& fvGeometry)
149
    {
150
151
152
153
        // for the o-scheme, the stencil is equal to the scv
        // index set of the dual grid's nodal index set
        stencil_ = &indexSet.nodalIndexSet().globalScvIndices();

Dennis Gläser's avatar
Dennis Gläser committed
154
        // number of interaction-volume-local (= node-local for o-scheme) scvs/scvf
155
156
157
158
        numFaces_ = indexSet.numFaces();
        const auto numLocalScvs = indexSet.numScvs();
        const auto numGlobalScvfs = indexSet.nodalIndexSet().numScvfs();

Dennis Gläser's avatar
Dennis Gläser committed
159
160
161
162
163
164
        // reserve memory for local entities
        elements_.clear();
        scvs_.clear();
        scvfs_.clear();
        localFaceData_.clear();
        dirichletData_.clear();
165
166
167
168
169
170
171
        elements_.reserve(numLocalScvs);
        scvs_.reserve(numLocalScvs);
        scvfs_.reserve(numFaces_);
        dirichletData_.reserve(numFaces_);
        localFaceData_.reserve(numGlobalScvfs);

        // set up quantities related to sub-control volumes
Dennis Gläser's avatar
Dennis Gläser committed
172
        const auto& scvIndices = indexSet.globalScvIndices();
173
174
        for (LocalIndexType scvIdxLocal = 0; scvIdxLocal < numLocalScvs; scvIdxLocal++)
        {
Dennis Gläser's avatar
Dennis Gläser committed
175
176
177
178
179
180
            scvs_.emplace_back(Helper(),
                               fvGeometry,
                               fvGeometry.scv( scvIndices[scvIdxLocal] ),
                               scvIdxLocal,
                               indexSet);
            elements_.emplace_back(fvGeometry.fvGridGeometry().element( scvIndices[scvIdxLocal] ));
181
182
183
184
        }

        // keep track of the number of unknowns etc
        numUnknowns_ = 0;
Dennis Gläser's avatar
Dennis Gläser committed
185
186
187
188
        numKnowns_ = numLocalScvs;

        // resize omega storage container
        wijk_.resize(numFaces_);
189
190
191
192

        // set up quantitites related to sub-control volume faces
        for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numFaces_; ++faceIdxLocal)
        {
Dennis Gläser's avatar
Dennis Gläser committed
193
194
            // get corresponding grid scvf
            const auto& scvf = fvGeometry.scvf(indexSet.scvfIdxGlobal(faceIdxLocal));
195

Dennis Gläser's avatar
Dennis Gläser committed
196
197
            // the neighboring scvs in local indices (order: 0 - inside scv, 1..n - outside scvs)
            const auto& neighborScvIndicesLocal = indexSet.neighboringLocalScvIndices(faceIdxLocal);
198
199

            // create local face data object for this face
Dennis Gläser's avatar
Dennis Gläser committed
200
            localFaceData_.emplace_back(faceIdxLocal, neighborScvIndicesLocal[0], scvf.index());
201
202
203
204

            // create iv-local scvf object
            if (scvf.boundary())
            {
Dennis Gläser's avatar
Dennis Gläser committed
205
                const auto bcTypes = problem.boundaryTypes(elements_[neighborScvIndicesLocal[0]], scvf);
206
207
208

                if (bcTypes.hasOnlyDirichlet())
                {
Dennis Gläser's avatar
Dennis Gläser committed
209
210
                    scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numKnowns_++, /*isDirichlet*/true);
                    dirichletData_.emplace_back(scvf.outsideScvIdx());
211
212
                }
                else
Dennis Gläser's avatar
Dennis Gläser committed
213
214
215
216
                    scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numUnknowns_++, /*isDirichlet*/false);

                // on boundary faces we will only need one inside omega
                wijk_[faceIdxLocal].resize(1);
217
218
219
            }
            else
            {
Dennis Gläser's avatar
Dennis Gläser committed
220
                scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numUnknowns_++, /*isDirichlet*/false);
221

Dennis Gläser's avatar
Dennis Gläser committed
222
223
224
                // we will need as many omegas as scvs around the face
                const auto numNeighborScvs = neighborScvIndicesLocal.size();
                wijk_[faceIdxLocal].resize(numNeighborScvs);
225

Dennis Gläser's avatar
Dennis Gläser committed
226
227
228
                // add local face data objects for the outside faces
                for (LocalIndexType i = 1; i < numNeighborScvs; ++i)
                {
229
                    // loop over scvfs in outside scv until we find the one coinciding with current scvf
Dennis Gläser's avatar
Dennis Gläser committed
230
                    const auto outsideLocalScvIdx = neighborScvIndicesLocal[i];
231
232
233
234
235
236
                    for (int coord = 0; coord < dim; ++coord)
                    {
                        if (indexSet.scvfIdxLocal(outsideLocalScvIdx, coord) == faceIdxLocal)
                        {
                            const auto globalScvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(outsideLocalScvIdx, coord);
                            const auto& flipScvf = fvGeometry.scvf(globalScvfIdx);
Dennis Gläser's avatar
Dennis Gläser committed
237
238
239
240
                            localFaceData_.emplace_back(faceIdxLocal,       // iv-local scvf idx
                                                        outsideLocalScvIdx, // iv-local scv index
                                                        i-1,                // scvf-local index in outside faces
                                                        flipScvf.index());  // global scvf index
241
242
243
244
245
246
                        }
                    }
                }
            }
        }

Dennis Gläser's avatar
Dennis Gläser committed
247
        // Resize local matrices
248
        A_.resize(numUnknowns_, numUnknowns_);
Dennis Gläser's avatar
Dennis Gläser committed
249
        B_.resize(numUnknowns_, numKnowns_);
250
        C_.resize(numFaces_, numUnknowns_);
251
252
    }

Dennis Gläser's avatar
Dennis Gläser committed
253
254
    //! returns the number of primary scvfs of this interaction volume
    std::size_t numFaces() const { return numFaces_; }
255

Dennis Gläser's avatar
Dennis Gläser committed
256
257
    //! returns the number of intermediate unknowns within this interaction volume
    std::size_t numUnknowns() const { return numUnknowns_; }
Dennis Gläser's avatar
Dennis Gläser committed
258

Dennis Gläser's avatar
Dennis Gläser committed
259
260
    //! returns the number of (in this context) known solution values within this interaction volume
    std::size_t numKnowns() const { return numKnowns_; }
Dennis Gläser's avatar
Dennis Gläser committed
261

Dennis Gläser's avatar
Dennis Gläser committed
262
263
    //! returns the number of scvs embedded in this interaction volume
    std::size_t numScvs() const { return scvs_.size(); }
264

Dennis Gläser's avatar
Dennis Gläser committed
265
266
    //! returns the number of scvfs embedded in this interaction volume
    std::size_t numScvfs() const { return scvfs_.size(); }
267

268
269
270
    //! returns the cell-stencil of this interaction volume
    const Stencil& stencil() const { return *stencil_; }

271
    //! returns the grid element corresponding to a given iv-local scv idx
272
    const Element& element(const LocalIndexType ivLocalScvIdx) const { return elements_[ivLocalScvIdx]; }
273

274
    //! returns the local scvf entity corresponding to a given iv-local scvf idx
275
    const LocalScvfType& localScvf(const LocalIndexType ivLocalScvfIdx) const { return scvfs_[ivLocalScvfIdx]; }
276

277
    //! returns the local scv entity corresponding to a given iv-local scv idx
278
    const LocalScvType& localScv(const LocalIndexType ivLocalScvfIdx) const { return scvs_[ivLocalScvfIdx]; }
279

280
    //! returns a reference to the container with the local face data
281
    const std::vector<LocalFaceData>& localFaceData() const { return localFaceData_; }
282

Dennis Gläser's avatar
Dennis Gläser committed
283
284
285
    //! returns a reference to the information container on Dirichlet BCs within this iv
    const std::vector<DirichletData>& dirichletData() const { return dirichletData_; }

286
    //! returns the matrix associated with face unknowns in local equation system
Dennis Gläser's avatar
Dennis Gläser committed
287
288
289
    const Matrix& A() const { return A_; }
    Matrix& A() { return A_; }

290
    //! returns the matrix associated with cell unknowns in local equation system
Dennis Gläser's avatar
Dennis Gläser committed
291
292
293
    const Matrix& B() const { return B_; }
    Matrix& B() { return B_; }

294
    //! returns the matrix associated with face unknowns in flux expressions
Dennis Gläser's avatar
Dennis Gläser committed
295
296
297
    const Matrix& C() const { return C_; }
    Matrix& C() { return C_; }

298
    //! returns container storing the transmissibilities for each face & coordinate
Dennis Gläser's avatar
Dennis Gläser committed
299
300
    const std::vector< std::vector< DimVector > >& omegas() const { return wijk_; }
    std::vector< std::vector< DimVector > >& omegas() { return wijk_; }
301

Dennis Gläser's avatar
Dennis Gläser committed
302
303
    //! returns the number of interaction volumes living around a vertex
    //! the mpfa-o scheme always constructs one iv per vertex
Dennis Gläser's avatar
Dennis Gläser committed
304
305
    template< class NodalIndexSet >
    static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet) { return 1; }
306

Dennis Gläser's avatar
Dennis Gläser committed
307
308
    //! adds the iv index sets living around a vertex to a given container
    //! and stores the the corresponding index in a map for each scvf
Dennis Gläser's avatar
Dennis Gläser committed
309
    template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet>
Dennis Gläser's avatar
Dennis Gläser committed
310
311
    static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
                                              ScvfIndexMap& scvfIndexMap,
Dennis Gläser's avatar
Dennis Gläser committed
312
                                              const NodalIndexSet& nodalIndexSet)
313
    {
Dennis Gläser's avatar
Dennis Gläser committed
314
315
        // the global index of the iv index set that is about to be created
        const auto curGlobalIndex = ivIndexSetContainer.size();
316

Dennis Gläser's avatar
Dennis Gläser committed
317
318
        // make the one index set for this node
        ivIndexSetContainer.emplace_back(nodalIndexSet);
319

Dennis Gläser's avatar
Dennis Gläser committed
320
321
322
        // store the index mapping
        for (const auto scvfIdx : nodalIndexSet.globalScvfIndices())
            scvfIndexMap[scvfIdx] = curGlobalIndex;
323
324
    }

325
private:
326
327
328
    // pointer to cell stencil (in iv index set)
    const Stencil* stencil_;

329
    // Variables defining the local scope
Dennis Gläser's avatar
Dennis Gläser committed
330
331
332
    std::vector<Element> elements_;
    std::vector<LocalScvType> scvs_;
    std::vector<LocalScvfType> scvfs_;
333
    std::vector<LocalFaceData> localFaceData_;
Dennis Gläser's avatar
Dennis Gläser committed
334
    std::vector<DirichletData> dirichletData_;
335

Dennis Gläser's avatar
Dennis Gläser committed
336
    // Matrices needed for computation of transmissibilities
Dennis Gläser's avatar
Dennis Gläser committed
337
338
339
    Matrix A_;
    Matrix B_;
    Matrix C_;
Dennis Gläser's avatar
Dennis Gläser committed
340
341
342
343
344
345
346
347

    // The omega factors are stored during assemble of local system
    std::vector< std::vector< DimVector > > wijk_;

    // sizes involved in the local system equations
    std::size_t numFaces_;
    std::size_t numUnknowns_;
    std::size_t numKnowns_;
348
349
350
351
352
};

} // end namespace

#endif