elementvolumevariables.hh 15.5 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 The local (stencil) volume variables class for cell centered mpfa models
23
24
25
26
 */
#ifndef DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH

27
28
29
#include <utility>
#include <dumux/common/properties.hh>

Dennis Gläser's avatar
Dennis Gläser committed
30
31
namespace Dumux
{
32
33

/*!
Dennis Gläser's avatar
Dennis Gläser committed
34
35
36
 * \ingroup CCMpfaDiscretization
 * \brief The local (stencil) volume variables class for cell centered mpfa models
 * \note The class is specilized for versions with and without caching
37
 */
38
template<class TypeTag, bool enableGridVolVarsCache>
Dennis Gläser's avatar
Dennis Gläser committed
39
class CCMpfaElementVolumeVariables;
40

Dennis Gläser's avatar
Dennis Gläser committed
41
42
43
44
45
/*!
 * \ingroup CCMpfaDiscretization
 * \brief The local (stencil) volume variables class for cell centered mpfa models with caching
 * \note the volume variables are stored for the whole grid view in the corresponding GridVolumeVariables class
 */
46
template<class TypeTag>
47
class CCMpfaElementVolumeVariables<TypeTag, /*enableGridVolVarsCache*/true>
48
49
50
{
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
Dennis Gläser's avatar
Dennis Gläser committed
51
52
    using Element = typename GridView::template Codim<0>::Entity;

53
54
55
    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
56
    using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
57
    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
Dennis Gläser's avatar
Dennis Gläser committed
58
    using GridIndexType = typename GridView::IndexSet::IndexType;
59
60
61

public:
    //! Constructor
62
63
    CCMpfaElementVolumeVariables(const GridVolumeVariables& gridVolVars)
    : gridVolVarsPtr_(&gridVolVars) {}
64

Dennis Gläser's avatar
Dennis Gläser committed
65
    //! operator for the access with an scv
66
    const VolumeVariables& operator [](const SubControlVolume& scv) const
67
    { return gridVolVars().volVars(scv.dofIndex()); }
68

Dennis Gläser's avatar
Dennis Gläser committed
69
70
    //! operator for the access with an index
    const VolumeVariables& operator [](const GridIndexType scvIdx) const
71
    { return gridVolVars().volVars(scvIdx); }
72

Dennis Gläser's avatar
Dennis Gläser committed
73
    //! precompute all volume variables in a stencil of an element - do nothing volVars: are cached
74
75
    void bind(const Element& element,
              const FVElementGeometry& fvGeometry,
Dennis Gläser's avatar
Dennis Gläser committed
76
77
              const SolutionVector& sol)
    {}
78

Dennis Gläser's avatar
Dennis Gläser committed
79
    //! precompute the volume variables of an element - do nothing: volVars are cached
80
81
    void bindElement(const Element& element,
                     const FVElementGeometry& fvGeometry,
Dennis Gläser's avatar
Dennis Gläser committed
82
83
                     const SolutionVector& sol)
    {}
84
85

    //! The global volume variables object we are a restriction of
86
87
    const GridVolumeVariables& gridVolVars() const
    { return *gridVolVarsPtr_; }
88
89

private:
90
    const GridVolumeVariables* gridVolVarsPtr_;
91
92
93
};


Dennis Gläser's avatar
Dennis Gläser committed
94
95
96
97
/*!
 * \ingroup CCMpfaDiscretization
 * \brief The local (stencil) volume variables class for cell centered tpfa models with caching
 */
98
template<class TypeTag>
99
class CCMpfaElementVolumeVariables<TypeTag, /*enableGridVolVarsCache*/false>
100
101
102
{
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
Dennis Gläser's avatar
Dennis Gläser committed
103
104
    using Element = typename GridView::template Codim<0>::Entity;

105
    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
106
    using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
107
    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
108
    using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
109
110
    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
Dennis Gläser's avatar
Dennis Gläser committed
111
    using GridIndexType = typename GridView::IndexSet::IndexType;
112
113
114
115


public:
    //! Constructor
116
117
    CCMpfaElementVolumeVariables(const GridVolumeVariables& gridVolVars)
    : gridVolVarsPtr_(&gridVolVars) {}
118

Dennis Gläser's avatar
Dennis Gläser committed
119
    //! Prepares the volume variables within the element stencil
120
121
122
123
    void bind(const Element& element,
              const FVElementGeometry& fvGeometry,
              const SolutionVector& sol)
    {
Dennis Gläser's avatar
Dennis Gläser committed
124
125
        clear();

126
        const auto& problem = gridVolVars().problem();
Dennis Gläser's avatar
Dennis Gläser committed
127
128
        const auto& fvGridGeometry = fvGeometry.fvGridGeometry();
        const auto& gridIvIndexSets = fvGridGeometry.gridInteractionVolumeIndexSets();
129
130

        // stencil information
Dennis Gläser's avatar
Dennis Gläser committed
131
132
133
        const auto globalI = fvGridGeometry.elementMapper().index(element);
        const auto& assemblyMapI = fvGridGeometry.connectivityMap()[globalI];
        const auto numVolVars = assemblyMapI.size() + 1;
134
135

        // resize local containers to the required size (for internal elements)
Dennis Gläser's avatar
Dennis Gläser committed
136
137
138
        volumeVariables_.resize(numVolVars);
        volVarIndices_.resize(numVolVars);
        unsigned int localIdx = 0;
139
140

        // update the volume variables of the element at hand
Dennis Gläser's avatar
Dennis Gläser committed
141
        const auto& scvI = fvGeometry.scv(globalI);
142
        volumeVariables_[localIdx].update(ElementSolution(sol[globalI]),
143
144
145
                                          problem,
                                          element,
                                          scvI);
146
        volVarIndices_[localIdx] = scvI.dofIndex();
147
148
149
        ++localIdx;

        // Update the volume variables of the neighboring elements
Dennis Gläser's avatar
Dennis Gläser committed
150
        for (auto&& dataJ : assemblyMapI)
151
        {
152
            const auto& elementJ = fvGridGeometry.element(dataJ.globalJ);
Dennis Gläser's avatar
Dennis Gläser committed
153
            const auto& scvJ = fvGeometry.scv(dataJ.globalJ);
154
            volumeVariables_[localIdx].update(ElementSolution(sol[dataJ.globalJ]),
155
156
157
                                              problem,
                                              elementJ,
                                              scvJ);
158
            volVarIndices_[localIdx] = scvJ.dofIndex();
159
160
161
            ++localIdx;
        }

Dennis Gläser's avatar
Dennis Gläser committed
162
        // maybe prepare boundary volume variables
Dennis Gläser's avatar
Dennis Gläser committed
163
164
        const auto maxNumBoundaryVolVars = maxNumBoundaryVolVars_(fvGeometry);
        if (maxNumBoundaryVolVars > 0)
Dennis Gläser's avatar
Dennis Gläser committed
165
        {
Dennis Gläser's avatar
Dennis Gläser committed
166
167
            volumeVariables_.reserve(numVolVars+maxNumBoundaryVolVars);
            volVarIndices_.reserve(numVolVars+maxNumBoundaryVolVars);
168
169

            // treat the BCs inside the element
Dennis Gläser's avatar
Dennis Gläser committed
170
            for (const auto& scvf : scvfs(fvGeometry))
171
            {
Dennis Gläser's avatar
Dennis Gläser committed
172
173
174
175
176
177
                // if we are not on a boundary, skip to the next scvf
                if (!scvf.boundary())
                    continue;

                const auto bcTypes = problem.boundaryTypes(element, scvf);

Dennis Gläser's avatar
Dennis Gläser committed
178
179
                // Only proceed on dirichlet boundaries. Fluxes across Neumann
                // boundaries are never computed - the user-defined flux is taken.
Dennis Gläser's avatar
Dennis Gläser committed
180
                if (bcTypes.hasOnlyDirichlet())
181
                {
Dennis Gläser's avatar
Dennis Gläser committed
182
183
                    // boundary volume variables
                    VolumeVariables dirichletVolVars;
184
                    dirichletVolVars.update(ElementSolution(problem.dirichlet(element, scvf)),
Dennis Gläser's avatar
Dennis Gläser committed
185
186
187
188
189
190
191
                                            problem,
                                            element,
                                            scvI);

                    volumeVariables_.emplace_back(std::move(dirichletVolVars));
                    volVarIndices_.push_back(scvf.outsideScvIdx());
                }
192
193
            }

194
            // Update boundary volume variables in the neighbors
Dennis Gläser's avatar
Dennis Gläser committed
195
            for (const auto& scvf : scvfs(fvGeometry))
196
            {
Dennis Gläser's avatar
Dennis Gläser committed
197
198
199
200
201
202
203
204
                if (fvGridGeometry.vertexUsesPrimaryInteractionVolume(scvf.vertexIndex()))
                {
                    const auto& nodalIndexSet = gridIvIndexSets.primaryIndexSet(scvf).nodalIndexSet();
                    // if present, insert boundary vol vars
                    if (nodalIndexSet.numBoundaryScvfs() > 0)
                        addBoundaryVolVars_(problem, fvGeometry, nodalIndexSet);
                }
                else
205
                {
Dennis Gläser's avatar
Dennis Gläser committed
206
207
208
209
                    const auto& nodalIndexSet = gridIvIndexSets.secondaryIndexSet(scvf).nodalIndexSet();
                    // if present, insert boundary vol vars
                    if (nodalIndexSet.numBoundaryScvfs() > 0)
                        addBoundaryVolVars_(problem, fvGeometry, nodalIndexSet);
210
                }
211
            }
212
        }
213

Dennis Gläser's avatar
Dennis Gläser committed
214
215
216
217
218
219
220
221
222
223
224
225
226
        // //! TODO Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
        // //! on additional DOFs not included in the discretization schemes' occupation pattern
        // const auto& additionalDofDependencies = problem.getAdditionalDofDependencies(globalI);
        // if (!additionalDofDependencies.empty())
        // {
        //     volumeVariables_.reserve(volumeVariables_.size() + additionalDofDependencies.size());
        //     volVarIndices_.reserve(volVarIndices_.size() + additionalDofDependencies.size());
        //     for (auto globalJ : additionalDofDependencies)
        //     {
        //         const auto& elementJ = fvGridGeometry.element(globalJ);
        //         const auto& scvJ = fvGeometry.scv(globalJ);

        //         VolumeVariables additionalVolVars;
227
        //         additionalVolVars.update(ElementSolution(elementJ, sol, fvGridGeometry),
Dennis Gläser's avatar
Dennis Gläser committed
228
229
230
231
232
233
234
235
        //                                  problem,
        //                                  elementJ,
        //                                  scvJ);

        //         volumeVariables_.emplace_back(std::move(additionalVolVars));
        //         volVarIndices_.push_back(globalJ);
        //     }
        // }
236
237
    }

Dennis Gläser's avatar
Dennis Gläser committed
238
    //! Prepares the volume variables of an element
239
240
241
242
    void bindElement(const Element& element,
                     const FVElementGeometry& fvGeometry,
                     const SolutionVector& sol)
    {
Dennis Gläser's avatar
Dennis Gläser committed
243
244
        clear();

Dennis Gläser's avatar
Dennis Gläser committed
245
        auto eIdx = fvGeometry.fvGridGeometry().elementMapper().index(element);
246
247
248
249
        volumeVariables_.resize(1);
        volVarIndices_.resize(1);

        // update the volume variables of the element
Dennis Gläser's avatar
Dennis Gläser committed
250
        const auto& scv = fvGeometry.scv(eIdx);
251
        volumeVariables_[0].update(ElementSolution(sol[scv.dofIndex()]),
252
                                   gridVolVars().problem(),
253
254
                                   element,
                                   scv);
255
        volVarIndices_[0] = scv.dofIndex();
256
257
    }

Dennis Gläser's avatar
Dennis Gläser committed
258
    //! access operator with scv
259
    const VolumeVariables& operator [](const SubControlVolume& scv) const
260
    { return volumeVariables_[getLocalIdx_(scv.dofIndex())]; }
261

Dennis Gläser's avatar
Dennis Gläser committed
262
    //! access operator with scv
263
    VolumeVariables& operator [](const SubControlVolume& scv)
264
    { return volumeVariables_[getLocalIdx_(scv.dofIndex())]; }
265

Dennis Gläser's avatar
Dennis Gläser committed
266
267
    //! access operator with scv index
    const VolumeVariables& operator [](GridIndexType scvIdx) const
268
269
    { return volumeVariables_[getLocalIdx_(scvIdx)]; }

Dennis Gläser's avatar
Dennis Gläser committed
270
271
    //! access operator with scv index
    VolumeVariables& operator [](GridIndexType scvIdx)
272
273
274
    { return volumeVariables_[getLocalIdx_(scvIdx)]; }

    //! The global volume variables object we are a restriction of
275
276
    const GridVolumeVariables& gridVolVars() const
    { return *gridVolVarsPtr_; }
277

Dennis Gläser's avatar
Dennis Gläser committed
278
279
280
281
282
283
284
    //! Clear all local storage
    void clear()
    {
        volVarIndices_.clear();
        volumeVariables_.clear();
    }

285
private:
286
    const GridVolumeVariables* gridVolVarsPtr_;
287

Dennis Gläser's avatar
Dennis Gläser committed
288
289
290
291
    // Computes how many boundary vol vars come into play for flux calculations
    // on this element. This number here is probably always higher than the actually
    // needed number of volume variables. However, memory is not an issue for the global
    // caching being deactivated and we want to make sure we reserve enough memory here.
Dennis Gläser's avatar
Dennis Gläser committed
292
293
    // TODO: What about non-symmetric schemes? Is there a better way for estimating this?
    std::size_t maxNumBoundaryVolVars_(const FVElementGeometry& fvGeometry)
294
    {
Dennis Gläser's avatar
Dennis Gläser committed
295
296
297
298
299
        const auto& fvGridGeometry = fvGeometry.fvGridGeometry();
        const auto& gridIvIndexSets = fvGridGeometry.gridInteractionVolumeIndexSets();

        std::size_t numBoundaryVolVars = 0;
        for (const auto& scvf : scvfs(fvGeometry))
300
        {
Dennis Gläser's avatar
Dennis Gläser committed
301
302
303
304
            if (fvGridGeometry.vertexUsesPrimaryInteractionVolume(scvf.vertexIndex()))
                numBoundaryVolVars += gridIvIndexSets.primaryIndexSet(scvf).nodalIndexSet().numBoundaryScvfs();
            else
                numBoundaryVolVars += gridIvIndexSets.secondaryIndexSet(scvf).nodalIndexSet().numBoundaryScvfs();
305
306
        }

Dennis Gläser's avatar
Dennis Gläser committed
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
        return numBoundaryVolVars;
    }

    //! adds the volume variables from a given nodal index set to the local containers
    template<class NodalIndexSet>
    void addBoundaryVolVars_(const Problem& problem, const FVElementGeometry& fvGeometry, const NodalIndexSet& nodalIndexSet)
    {
        // check each scvf in the index set for boundary presence
        for (auto scvfIdx : nodalIndexSet.globalScvfIndices())
        {
            const auto& ivScvf = fvGeometry.scvf(scvfIdx);

            // only proceed for scvfs on the boundary and not in the inside element
            if (!ivScvf.boundary() || ivScvf.insideScvIdx() == volVarIndices_[0])
                continue;

            const auto insideScvIdx = ivScvf.insideScvIdx();
            const auto insideElement = fvGeometry.fvGridGeometry().element(insideScvIdx);
            const auto bcTypes = problem.boundaryTypes(insideElement, ivScvf);

Dennis Gläser's avatar
Dennis Gläser committed
327
328
            // Only proceed on dirichlet boundaries. Fluxes across Neumann
            // boundaries are never computed - the user-defined flux is taken.
Dennis Gläser's avatar
Dennis Gläser committed
329
330
331
332
            if (bcTypes.hasOnlyDirichlet())
            {
                VolumeVariables dirichletVolVars;
                const auto& ivScv = fvGeometry.scv(insideScvIdx);
333
                dirichletVolVars.update(ElementSolution(problem.dirichlet(insideElement, ivScvf)),
Dennis Gläser's avatar
Dennis Gläser committed
334
335
336
337
338
339
340
341
                                        problem,
                                        insideElement,
                                        ivScv);

                volumeVariables_.emplace_back(std::move(dirichletVolVars));
                volVarIndices_.push_back(ivScvf.outsideScvIdx());
            }
        }
342
343
    }

Dennis Gläser's avatar
Dennis Gläser committed
344
    //! map a global scv index to the local storage index
345
    int getLocalIdx_(const int volVarIdx) const
346
347
348
349
350
351
    {
        auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
        assert(it != volVarIndices_.end() && "Could not find the current volume variables for volVarIdx!");
        return std::distance(volVarIndices_.begin(), it);
    }

Dennis Gläser's avatar
Dennis Gläser committed
352
    std::vector<GridIndexType> volVarIndices_;
353
354
355
356
357
358
    std::vector<VolumeVariables> volumeVariables_;
};

} // end namespace

#endif