interactionvolume.hh 15.2 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
22
 * \ingroup CCMpfaDiscretization
 * \brief Classes 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
86
87
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 >;
    //! export the data handle type for this iv
    using DataHandle = InteractionVolumeDataHandle< TypeTag, InteractionVolumeType >;
88
};
89

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

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

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

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

Dennis Gläser's avatar
Dennis Gläser committed
116
117
118
    using IndexSet = typename TraitsType::IndexSet;
    using GridIndexType = typename TraitsType::GridIndexType;
    using LocalIndexType = typename TraitsType::LocalIndexType;
119

Dennis Gläser's avatar
Dennis Gläser committed
120
121
122
123
124
125
126
127
128
129
130
131
132
    //! 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
133

Dennis Gläser's avatar
Dennis Gläser committed
134
public:
135
136
137
    //! 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
138
139
140
141
    //! Types required in the assembly of the local eq system
    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
142

Dennis Gläser's avatar
Dennis Gläser committed
143
    //! Sets up the local scope for a given iv index set
144
145
146
    void setUpLocalScope(const IndexSet& indexSet,
                         const Problem& problem,
                         const FVElementGeometry& fvGeometry)
147
    {
Dennis Gläser's avatar
Dennis Gläser committed
148
        // number of interaction-volume-local (= node-local for o-scheme) scvs/scvf
149
150
151
152
        numFaces_ = indexSet.numFaces();
        const auto numLocalScvs = indexSet.numScvs();
        const auto numGlobalScvfs = indexSet.nodalIndexSet().numScvfs();

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

        // keep track of the number of unknowns etc
        numUnknowns_ = 0;
Dennis Gläser's avatar
Dennis Gläser committed
179
180
181
182
        numKnowns_ = numLocalScvs;

        // resize omega storage container
        wijk_.resize(numFaces_);
183
184
185
186

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

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

            // create local face data object for this face
Dennis Gläser's avatar
Dennis Gläser committed
194
            localFaceData_.emplace_back(faceIdxLocal, neighborScvIndicesLocal[0], scvf.index());
195
196
197
198

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

                if (bcTypes.hasOnlyDirichlet())
                {
Dennis Gläser's avatar
Dennis Gläser committed
203
204
                    scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numKnowns_++, /*isDirichlet*/true);
                    dirichletData_.emplace_back(scvf.outsideScvIdx());
205
206
                }
                else
Dennis Gläser's avatar
Dennis Gläser committed
207
208
209
210
                    scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numUnknowns_++, /*isDirichlet*/false);

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

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

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

Dennis Gläser's avatar
Dennis Gläser committed
241
        // Resize local matrices
242
        A_.resize(numUnknowns_, numUnknowns_);
Dennis Gläser's avatar
Dennis Gläser committed
243
        B_.resize(numUnknowns_, numKnowns_);
244
        C_.resize(numFaces_, numUnknowns_);
245
246
    }

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

Dennis Gläser's avatar
Dennis Gläser committed
250
251
    //! 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
252

Dennis Gläser's avatar
Dennis Gläser committed
253
254
    //! 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
255

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

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

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

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

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

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

Dennis Gläser's avatar
Dennis Gläser committed
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
    //! returns a reference to the information container on Dirichlet BCs within this iv
    const std::vector<DirichletData>& dirichletData() const { return dirichletData_; }

    //! return functions for matrices involved in local system of equations
    const Matrix& A() const { return A_; }
    Matrix& A() { return A_; }

    const Matrix& B() const { return B_; }
    Matrix& B() { return B_; }

    const Matrix& C() const { return C_; }
    Matrix& C() { return C_; }

    const std::vector< std::vector< DimVector > >& omegas() const { return wijk_; }
    std::vector< std::vector< DimVector > >& omegas() { return wijk_; }
289

Dennis Gläser's avatar
Dennis Gläser committed
290
291
    //! 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
292
293
    template< class NodalIndexSet >
    static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet) { return 1; }
294

Dennis Gläser's avatar
Dennis Gläser committed
295
296
    //! 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
297
    template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet>
Dennis Gläser's avatar
Dennis Gläser committed
298
299
    static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
                                              ScvfIndexMap& scvfIndexMap,
Dennis Gläser's avatar
Dennis Gläser committed
300
                                              const NodalIndexSet& nodalIndexSet)
301
    {
Dennis Gläser's avatar
Dennis Gläser committed
302
303
        // the global index of the iv index set that is about to be created
        const auto curGlobalIndex = ivIndexSetContainer.size();
304

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

Dennis Gläser's avatar
Dennis Gläser committed
308
309
310
        // store the index mapping
        for (const auto scvfIdx : nodalIndexSet.globalScvfIndices())
            scvfIndexMap[scvfIdx] = curGlobalIndex;
311
312
    }

313
private:
314
    // Variables defining the local scope
Dennis Gläser's avatar
Dennis Gläser committed
315
316
317
    std::vector<Element> elements_;
    std::vector<LocalScvType> scvs_;
    std::vector<LocalScvfType> scvfs_;
318
    std::vector<LocalFaceData> localFaceData_;
Dennis Gläser's avatar
Dennis Gläser committed
319
    std::vector<DirichletData> dirichletData_;
320

Dennis Gläser's avatar
Dennis Gläser committed
321
    // Matrices needed for computation of transmissibilities
Dennis Gläser's avatar
Dennis Gläser committed
322
323
324
    Matrix A_;
    Matrix B_;
    Matrix C_;
Dennis Gläser's avatar
Dennis Gläser committed
325
326
327
328
329
330
331
332

    // 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_;
333
334
335
336
337
};

} // end namespace

#endif