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>
36
#include <dumux/discretization/cellcentered/mpfa/interactionvolume.hh>
Dennis Gläser's avatar
Dennis Gläser committed
37
#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
Dennis Gläser's avatar
Dennis Gläser committed
38
39
#include <dumux/discretization/cellcentered/mpfa/localfacedata.hh>
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
40
41

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

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

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

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

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

Dennis Gläser's avatar
Dennis Gläser committed
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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 >;
87
88
    //! export the type used for the iv-stencils
    using Stencil = std::vector< GridIndexType >;
Dennis Gläser's avatar
Dennis Gläser committed
89
90
    //! export the data handle type for this iv
    using DataHandle = InteractionVolumeDataHandle< TypeTag, InteractionVolumeType >;
91
};
92

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

105
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
Dennis Gläser's avatar
Dennis Gläser committed
106
    using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
107
108
    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
109
110

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

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

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

Dennis Gläser's avatar
Dennis Gläser committed
124
125
126
127
128
129
130
131
132
133
134
135
136
    //! 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
137

Dennis Gläser's avatar
Dennis Gläser committed
138
public:
139
140
141
    //! 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
142
143
144
    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
145

Dennis Gläser's avatar
Dennis Gläser committed
146
    //! Sets up the local scope for a given iv index set
147
148
149
    void setUpLocalScope(const IndexSet& indexSet,
                         const Problem& problem,
                         const FVElementGeometry& fvGeometry)
150
    {
151
152
153
154
        // 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
155
        // number of interaction-volume-local (= node-local for o-scheme) scvs/scvf
156
157
158
159
        numFaces_ = indexSet.numFaces();
        const auto numLocalScvs = indexSet.numScvs();
        const auto numGlobalScvfs = indexSet.nodalIndexSet().numScvfs();

Dennis Gläser's avatar
Dennis Gläser committed
160
161
162
163
164
165
        // reserve memory for local entities
        elements_.clear();
        scvs_.clear();
        scvfs_.clear();
        localFaceData_.clear();
        dirichletData_.clear();
166
167
168
169
170
171
172
        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
173
        const auto& scvIndices = indexSet.globalScvIndices();
174
175
        for (LocalIndexType scvIdxLocal = 0; scvIdxLocal < numLocalScvs; scvIdxLocal++)
        {
Dennis Gläser's avatar
Dennis Gläser committed
176
177
178
179
180
181
            scvs_.emplace_back(Helper(),
                               fvGeometry,
                               fvGeometry.scv( scvIndices[scvIdxLocal] ),
                               scvIdxLocal,
                               indexSet);
            elements_.emplace_back(fvGeometry.fvGridGeometry().element( scvIndices[scvIdxLocal] ));
182
183
184
185
        }

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

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

        // 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
194
195
            // get corresponding grid scvf
            const auto& scvf = fvGeometry.scvf(indexSet.scvfIdxGlobal(faceIdxLocal));
196

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

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

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

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

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

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

Dennis Gläser's avatar
Dennis Gläser committed
227
228
229
                // add local face data objects for the outside faces
                for (LocalIndexType i = 1; i < numNeighborScvs; ++i)
                {
230
                    // loop over scvfs in outside scv until we find the one coinciding with current scvf
Dennis Gläser's avatar
Dennis Gläser committed
231
                    const auto outsideLocalScvIdx = neighborScvIndicesLocal[i];
232
233
234
235
236
237
                    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
238
239
240
241
                            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
242
243
244
245
246
247
                        }
                    }
                }
            }
        }

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

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

Dennis Gläser's avatar
Dennis Gläser committed
257
258
    //! 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
259

Dennis Gläser's avatar
Dennis Gläser committed
260
261
    //! 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
262

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

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

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

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

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

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

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

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

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

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

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

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

Dennis Gläser's avatar
Dennis Gläser committed
303
304
    //! 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
305
306
    template< class NodalIndexSet >
    static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet) { return 1; }
307

Dennis Gläser's avatar
Dennis Gläser committed
308
309
    //! 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
310
    template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet>
Dennis Gläser's avatar
Dennis Gläser committed
311
312
    static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
                                              ScvfIndexMap& scvfIndexMap,
Dennis Gläser's avatar
Dennis Gläser committed
313
                                              const NodalIndexSet& nodalIndexSet)
314
    {
Dennis Gläser's avatar
Dennis Gläser committed
315
316
        // the global index of the iv index set that is about to be created
        const auto curGlobalIndex = ivIndexSetContainer.size();
317

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

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

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

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

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

    // 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_;
349
350
351
352
353
};

} // end namespace

#endif