elementvolumevariables.hh 11.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
// -*- 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
 * \brief The local (stencil) volume variables class for cell centered models
 */
#ifndef DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH

#include <dumux/implicit/properties.hh>

namespace Dumux
{

/*!
 * \ingroup ImplicitModel
 * \brief Base class for the volume variables vector
 */
template<class TypeTag, bool enableGlobalVolVarsCache>
class CCMpfaElementVolumeVariables
{};

// specialization in case of storing the volume variables globally
template<class TypeTag>
class CCMpfaElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
{
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
    using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
    using IndexType = typename GridView::IndexSet::IndexType;

    static const int dim = GridView::dimension;
    using Element = typename GridView::template Codim<0>::Entity;

public:
    //! Constructor
    CCMpfaElementVolumeVariables(const GlobalVolumeVariables& globalVolVars)
    : globalVolVarsPtr_(&globalVolVars) {}

    const VolumeVariables& operator [](const SubControlVolume& scv) const
    { return globalVolVars().volVars(scv.index()); }

    // operator for the access with an index
    // needed for cc methods for the access to the boundary volume variables
    const VolumeVariables& operator [](const IndexType scvIdx) const
    { return globalVolVars().volVars(scvIdx); }

    // For compatibility reasons with the case of not storing the vol vars.
    // function to be called before assembling an element, preparing the vol vars within the stencil
    void bind(const Element& element,
              const FVElementGeometry& fvGeometry,
72
              const SolutionVector& sol) {}
73
74
75
76

    // function to prepare the vol vars within the element
    void bindElement(const Element& element,
                     const FVElementGeometry& fvGeometry,
77
                     const SolutionVector& sol) {}
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

    //! The global volume variables object we are a restriction of
    const GlobalVolumeVariables& globalVolVars() const
    { return *globalVolVarsPtr_; }

private:
    const GlobalVolumeVariables* globalVolVarsPtr_;
};


// Specialization when the current volume variables are not stored
template<class TypeTag>
class CCMpfaElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false>
{
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
    using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using IndexType = typename GridView::IndexSet::IndexType;

    static const int dim = GridView::dimension;
    using Element = typename GridView::template Codim<0>::Entity;

public:

    //! Constructor
    CCMpfaElementVolumeVariables(const GlobalVolumeVariables& globalVolVars)
    : globalVolVarsPtr_(&globalVolVars) {}

    // Binding of an element, prepares the volume variables within the element stencil
    // called by the local jacobian to prepare element assembly
    void bind(const Element& element,
              const FVElementGeometry& fvGeometry,
              const SolutionVector& sol)
    {
116
117
        const auto& problem = globalVolVars().problem_();
        auto eIdx = problem.elementMapper().index(element);
118
119

        // stencil information
120
        const auto& neighborStencil = problem.model().stencils(element).neighborStencil();
121
122
123
124
125
126
127
128
129
        const auto numDofs = neighborStencil.size() + 1;

        // resize local containers to the required size (for internal elements)
        volumeVariables_.resize(numDofs);
        volVarIndices_.resize(numDofs);
        int localIdx = 0;

        // update the volume variables of the element at hand
        auto&& scvI = fvGeometry.scv(eIdx);
130
        volumeVariables_[localIdx].update(sol[eIdx], problem, element, scvI);
131
132
133
134
135
136
137
138
        volVarIndices_[localIdx] = scvI.index();
        ++localIdx;

        // Update the volume variables of the neighboring elements
        for (auto globalJ : neighborStencil)
        {
            const auto& elementJ = fvGeometry.globalFvGeometry().element(globalJ);
            auto&& scvJ = fvGeometry.scv(globalJ);
139
            volumeVariables_[localIdx].update(sol[globalJ], problem, elementJ, scvJ);
140
141
142
143
            volVarIndices_[localIdx] = scvJ.index();
            ++localIdx;
        }

144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
        std::size_t neighborScvfEstimate = neighborStencil.size()*dim*(2*dim-2);
        std::vector<IndexType> finishedBoundaries;
        finishedBoundaries.reserve(neighborScvfEstimate);

        // if the element is connected to a boundary, additionally treat the BC
        if (element.hasBoundaryIntersections())
        {
            // reserve enough space for the boundary volume variables
            volumeVariables_.reserve(numDofs+neighborScvfEstimate);
            volVarIndices_.reserve(numDofs+neighborScvfEstimate);

            for (auto&& scvf : scvfs(fvGeometry))
            {
                // if we are not on a boundary, skip to the next scvf
                if (!scvf.boundary())
                    continue;

                // boundary volume variables
                VolumeVariables dirichletVolVars;
                const auto dirichletPriVars = problem.dirichlet(element, scvf);
                dirichletVolVars.update(dirichletPriVars, problem, element, scvI);
                volumeVariables_.emplace_back(std::move(dirichletVolVars));

                // boundary vol var index
                auto bVolVarIdx = scvf.outsideScvIdx();
                volVarIndices_.push_back(bVolVarIdx);
                finishedBoundaries.push_back(bVolVarIdx);
            }
        }

        // Update boundary volume variables in the neighbors
        const auto& globalFvGeometry = problem.model().globalFvGeometry();
176
177
        for (auto&& scvf : scvfs(fvGeometry))
        {
178
179
180
181
            // reserve enough space for the boundary volume variables
            volumeVariables_.reserve(numDofs+neighborScvfEstimate);
            volVarIndices_.reserve(numDofs+neighborScvfEstimate);

182
            // skip the rest if the scvf does not touch a boundary
183
            const bool boundary = globalFvGeometry.scvfTouchesBoundary(scvf);
184
            if (!boundary)
185
186
                continue;

187
            // loop over all the scvfs in the interaction region
188
            const auto& ivSeed = globalFvGeometry.boundaryInteractionVolumeSeed(scvf);
189
190
191
192
            for (auto scvfIdx : ivSeed.globalScvfIndices())
            {
                auto&& ivScvf = fvGeometry.scvf(scvfIdx);

193
                if (!ivScvf.boundary() || contains(ivScvf.outsideScvIdx(), finishedBoundaries))
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
                    continue;

                // that means we are on a not yet handled boundary scvf
                auto insideScvIdx = ivScvf.insideScvIdx();
                auto&& ivScv = fvGeometry.scv(insideScvIdx);
                auto ivElement = globalFvGeometry.element(insideScvIdx);

                // boundary volume variables
                VolumeVariables dirichletVolVars;
                const auto dirichletPriVars = problem.dirichlet(ivElement, ivScvf);
                dirichletVolVars.update(dirichletPriVars, problem, ivElement, ivScv);
                volumeVariables_.emplace_back(std::move(dirichletVolVars));

                // boundary vol var index
                auto bVolVarIdx = ivScvf.outsideScvIdx();
                volVarIndices_.push_back(bVolVarIdx);
                finishedBoundaries.push_back(bVolVarIdx);
            }
212
        }
213
214
215
216

        // free unused memory
        volumeVariables_.shrink_to_fit();
        volVarIndices_.shrink_to_fit();
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
    }

    // Binding of an element, prepares only the volume variables of the element
    // specialization for cc models
    void bindElement(const Element& element,
                     const FVElementGeometry& fvGeometry,
                     const SolutionVector& sol)
    {
        auto eIdx = globalVolVars().problem_().elementMapper().index(element);
        volumeVariables_.resize(1);
        volVarIndices_.resize(1);

        // update the volume variables of the element
        auto&& scv = fvGeometry.scv(eIdx);
        volumeVariables_[0].update(sol[eIdx], globalVolVars().problem_(), element, scv);
        volVarIndices_[0] = scv.index();
    }

    const VolumeVariables& operator [](const SubControlVolume& scv) const
    { return volumeVariables_[getLocalIdx_(scv.index())]; }

    VolumeVariables& operator [](const SubControlVolume& scv)
    { return volumeVariables_[getLocalIdx_(scv.index())]; }

    const VolumeVariables& operator [](IndexType scvIdx) const
    { return volumeVariables_[getLocalIdx_(scvIdx)]; }

    VolumeVariables& operator [](IndexType scvIdx)
    { return volumeVariables_[getLocalIdx_(scvIdx)]; }

    //! The global volume variables object we are a restriction of
    const GlobalVolumeVariables& globalVolVars() const
    { return *globalVolVarsPtr_; }

private:
    const GlobalVolumeVariables* globalVolVarsPtr_;

    const int getLocalIdx_(const int volVarIdx) const
    {
        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);
    }

261
262
263
    const bool contains(const IndexType idx,
                        const std::vector<IndexType>& indices) const
    { return std::find(indices.begin(), indices.end(), idx) != indices.end(); }
264

265
266
267
268
269
270
271
    std::vector<IndexType> volVarIndices_;
    std::vector<VolumeVariables> volumeVariables_;
};

} // end namespace

#endif