interactionvolume.hh 29.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
// -*- 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 Base classes for interaction volume of mpfa methods.
 */
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH

#include <dumux/common/math.hh>
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
28
#include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh>
29
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
Dennis Gläser's avatar
Dennis Gläser committed
30
#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
31
32

#include "localsubcontrolentities.hh"
Dennis Gläser's avatar
Dennis Gläser committed
33
#include "interactionvolumeindexset.hh"
34
35
36

namespace Dumux
{
37
//! Specialization of the interaction volume traits class for the mpfa-o method
38
template<class TypeTag>
39
class CCMpfaOInteractionVolumeTraits : public CCMpfaInteractionVolumeTraitsBase<TypeTag>
40
{
41
    using BaseTraits = CCMpfaInteractionVolumeTraitsBase<TypeTag>;
Dennis Gläser's avatar
Dennis Gläser committed
42
    using NodalIndexSet = typename CCMpfaDualGridIndexSet<TypeTag>::NodalIndexSet;
43

44
public:
Dennis Gläser's avatar
Dennis Gläser committed
45
    using SecondaryInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>;
46

Dennis Gläser's avatar
Dennis Gläser committed
47
48
49
    using typename BaseTraits::DynamicLocalIndexContainer;
    using typename BaseTraits::DynamicGlobalIndexContainer;
    using IndexSet = CCMpfaOInteractionVolumeIndexSet<NodalIndexSet, DynamicGlobalIndexContainer, DynamicLocalIndexContainer>;
50

Dennis Gläser's avatar
Dennis Gläser committed
51
52
53
54
55
56
57
58
59
    // The matrix & vector types used in the interaction
    // volume are actually the dynamic types in the o-scheme
    using typename BaseTraits::DynamicVector;
    using typename BaseTraits::DynamicMatrix;
    using Vector = DynamicVector;
    using Matrix = DynamicMatrix;

    using LocalScvType = CCMpfaOInteractionVolumeLocalScv<TypeTag, IndexSet>;
    using LocalScvfType = CCMpfaOInteractionVolumeLocalScvf<TypeTag>;
60
};
61
62

//! Forward declaration of the mpfa-o interaction volume
Dennis Gläser's avatar
Dennis Gläser committed
63
template<class TypeTag, class Traits>
64
class CCMpfaOInteractionVolume;
65

66
67
/*!
 * \ingroup Mpfa
68
69
70
71
 * \brief Base class for the interaction volumes of the mpfa-o method.
 *        We introduce one more level of inheritance here because the o-method with
 *        full pressure support uses the mpfa-o interaction volume but uses a different
 *        traits class.
72
73
 */
template<class TypeTag>
74
75
class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>
          : public CCMpfaOInteractionVolume<TypeTag,
Dennis Gläser's avatar
Dennis Gläser committed
76
                                            CCMpfaOInteractionVolumeTraits<TypeTag>>
77
{
78
    using TraitsType = CCMpfaOInteractionVolumeTraits<TypeTag>;
79
public:
80
81
    // state the traits class type
    using Traits = TraitsType;
82
83
};

Dennis Gläser's avatar
Dennis Gläser committed
84
template<class TypeTag, class Traits>
85
class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase<TypeTag, Traits>
86
{
87
    using ParentType = CCMpfaInteractionVolumeBase<TypeTag, Traits>;
88
89
90
91
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
92
    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
93
94
95
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
Dennis Gläser's avatar
Dennis Gläser committed
96
    using DualGridNodalIndexSet = typename CCMpfaDualGridIndexSet<TypeTag>::NodalIndexSet;
97

98
    static const int dim = GridView::dimension;
99
    static const int dimWorld = GridView::dimensionworld;
100
    using Element = typename GridView::template Codim<0>::Entity;
101
102
    using DimVector = Dune::FieldVector<Scalar, dim>;
    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
103

Dennis Gläser's avatar
Dennis Gläser committed
104
105
    using Vector = typename Traits::DynamicVector;
    using Matrix = typename Traits::DynamicMatrix;
106
    using Tensor = typename Traits::Tensor;
107

108
109
    using LocalScvType = typename Traits::LocalScvType;
    using LocalScvfType = typename Traits::LocalScvfType;
110

Dennis Gläser's avatar
Dennis Gläser committed
111
112
113
114
115
    using IndexSet = typename Traits::IndexSet;
    using LocalIndexContainer = typename Traits::DynamicLocalIndexContainer;
    using LocalIndexType = typename LocalIndexContainer::value_type;
    using GlobalIndexContainer = typename Traits::DynamicGlobalIndexContainer;
    using DataHandle = typename Traits::DataHandle;
116
public:
Dennis Gläser's avatar
Dennis Gläser committed
117

118
    using typename ParentType::GlobalLocalFaceDataPair;
119
    using typename ParentType::LocalFaceData;
Dennis Gläser's avatar
Dennis Gläser committed
120
121
122
123
124
125
126
127
128
    using typename ParentType::DirichletDataContainer;

    //! Sets up the local scope for a given seed
    //! This function has to be called before using the IV!
    void bind(const IndexSet& indexSet,
              const Problem& problem,
              const FVElementGeometry& fvGeometry,
              const ElementVolumeVariables& elemVolVars,
              DataHandle& dataHandle)
129
    {
Dennis Gläser's avatar
Dennis Gläser committed
130
131
132
133
        problemPtr_ = &problem;
        fvGeometryPtr_ = &fvGeometry;
        elemVolVarsPtr_ = &elemVolVars;
        indexSetPtr_ = &indexSet;
134

Dennis Gläser's avatar
Dennis Gläser committed
135
136
        // set up the local scope of this interaction volume
        setLocalScope_(dataHandle);
137
138
    }

139
140
141
142
143
144
145
146
147
148
149
150
151
    //! Sets only the pointers to the local views
    //! Using the IV afterwards requires having called bind once before!
    //! Calling this with an fvGeometry or elemVolVars that differ from the
    //! ones that have been passed when calling bind leads to undefined behaviour
    void resetPointers(const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars)
    {
        fvGeometryPtr_ = &fvGeometry;
        elemVolVarsPtr_ = &elemVolVars;

        for (unsigned i = 0; i < globalLocalScvfPairedData_.size(); ++i)
            globalLocalScvfPairedData_[i].first = &fvGeometry.scvf(globalScvfIndices_[i]);
    }

Dennis Gläser's avatar
Dennis Gläser committed
152
    //! solves for the transmissibilities subject to a given tensor
153
    template<typename GetTensorFunction>
Dennis Gläser's avatar
Dennis Gläser committed
154
    void solveLocalSystem(const GetTensorFunction& getTensor, DataHandle& dataHandle)
155
156
    {
        // if only dirichlet faces are present, assemble T_ directly
Dennis Gläser's avatar
Dennis Gläser committed
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
        if (numUnknowns_ == 0)
            assemblePureDirichletSystem_(getTensor, dataHandle.T());
        else
        {
            // assemble
            assembleLocalMatrices_(getTensor);

            // solve
            A_.invert();

            // T = C*A^-1*B + D
            auto& T = dataHandle.T();
            T = multiplyMatrices(C_.rightmultiply(A_), B_);
            T += D_;

            // store A-1B only when gradient reconstruction is necessary
            if (requireABMatrix_())
                dataHandle.AB() = B_.leftmultiply(A_);
        }

        // // set vol vars stencil & positions pointer in handle
178
        dataHandle.setVolVarsStencilPointer(indexSet().globalScvIndices());
Dennis Gläser's avatar
Dennis Gläser committed
179
180
181
182
183
        dataHandle.setDirichletDataPointer(dirichletData_);

        // on surface grids, additionally prepare the outside transmissibilities
        if (dim < dimWorld)
            computeOutsideTransmissibilities_(dataHandle);
184
185
    }

186
    //! Gets the transmissibilities for a sub-control volume face within the interaction volume.
187
188
    //! specialization for dim == dimWorld
    template<int d = dim, int dw = dimWorld>
Dennis Gläser's avatar
Dennis Gläser committed
189
190
191
    typename std::enable_if< (d == dw), const Vector& >::type
    getTransmissibilities(const SubControlVolumeFace& scvf, const LocalFaceData& localFaceData, const DataHandle& dataHandle) const
    { return dataHandle.T()[localFaceData.localScvfIndex]; }
192

193
194
    //! Gets the transmissibilities for a sub-control volume face within the interaction volume.
    //! specialization for dim < dimWorld.
195
    template<int d = dim, int dw = dimWorld>
Dennis Gläser's avatar
Dennis Gläser committed
196
197
    typename std::enable_if< (d < dw), const Vector& >::type
    getTransmissibilities(const SubControlVolumeFace& scvf, const LocalFaceData& localFaceData, const DataHandle& dataHandle) const
198
    {
199
        // If we come from the inside, simply return tij
200
        if (!localFaceData.isOutside)
Dennis Gläser's avatar
Dennis Gläser committed
201
202
203
204
            return dataHandle.T()[localFaceData.localScvfIndex];
        else
            return dataHandle.outsideTij()[this->findIndexInVector(outsideScvfIndices_, scvf.index())];
    }
205

Dennis Gläser's avatar
Dennis Gläser committed
206
207
208
    //! obtain the local data object for a given global scvf
    const LocalFaceData& getLocalFaceData(const SubControlVolumeFace& scvf) const
    { return globalLocalScvfPairedData_[this->findIndexInVector(globalScvfIndices_, scvf.index())].second; }
209

Dennis Gläser's avatar
Dennis Gläser committed
210
211
212
    //! returns the vector of data pairs storing the local face data for each global scvf
    const std::vector<GlobalLocalFaceDataPair>& globalLocalScvfPairedData() const
    { return globalLocalScvfPairedData_; }
213

Dennis Gläser's avatar
Dennis Gläser committed
214
215
216
    //! returns the grid element corresponding to a given local (in the iv) scv idx
    const Element& element(const LocalIndexType localScvIdx) const
    { return elements_[localScvIdx]; }
217

Dennis Gläser's avatar
Dennis Gläser committed
218
219
220
    //! returns the local scvf entity corresponding to a given local (in the iv) scvf idx
    const LocalScvfType& localScvf(const LocalIndexType localScvfIdx) const
    { return scvfs_[localScvfIdx]; }
221

Dennis Gläser's avatar
Dennis Gläser committed
222
223
224
    //! returns the local scv entity corresponding to a given local (in the iv) scv idx
    const LocalScvType& localScv(const LocalIndexType localScvIdx) const
    { return scvs_[localScvIdx]; }
225

Dennis Gläser's avatar
Dennis Gläser committed
226
227
228
    //! returns a reference to the problem to be solved
    const Problem& problem() const
    { return *problemPtr_; }
229

Dennis Gläser's avatar
Dennis Gläser committed
230
231
232
    //! returns a reference to the fvGeometry object
    const FVElementGeometry& fvGeometry() const
    { return *fvGeometryPtr_; }
233

Dennis Gläser's avatar
Dennis Gläser committed
234
235
236
    //! returns a reference to the element volume variables
    const ElementVolumeVariables& elemVolVars() const
    { return *elemVolVarsPtr_; }
237

Dennis Gläser's avatar
Dennis Gläser committed
238
239
240
    //! returns a reference to the corresponding iv index set
    const IndexSet& indexSet() const
    { return *indexSetPtr_; }
241

Dennis Gläser's avatar
Dennis Gläser committed
242
    const DirichletDataContainer& dirichletData() const { return dirichletData_; }
243

Dennis Gläser's avatar
Dennis Gläser committed
244
245
246
247
    //! returns the number of interaction volumes living around a vertex
    //! the mpfa-o scheme always constructs one iv per vertex
    static std::size_t numInteractionVolumesAtVertex(const DualGridNodalIndexSet& nodalIndexSet)
    { return 1; }
248

Dennis Gläser's avatar
Dennis Gläser committed
249
250
251
252
253
254
    //! 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
    template<class IvIndexSetContainer, class ScvfIndexMap>
    static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
                                              ScvfIndexMap& scvfIndexMap,
                                              const DualGridNodalIndexSet& nodalIndexSet)
255
    {
Dennis Gläser's avatar
Dennis Gläser committed
256
257
        // the global index of the iv index set that is about to be created
        const auto curGlobalIndex = ivIndexSetContainer.size();
258

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

Dennis Gläser's avatar
Dennis Gläser committed
262
263
264
        // store the index mapping
        for (const auto scvfIdx : nodalIndexSet.globalScvfIndices())
            scvfIndexMap[scvfIdx] = curGlobalIndex;
265
266
    }

267
private:
Dennis Gläser's avatar
Dennis Gläser committed
268
269
270
271
272
273
    //! returns a boolean whether or not the AB matrix has to be passed to the handles
    bool requireABMatrix_() const
    {
        static const bool requireAB = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity) || dim < dimWorld;
        return requireAB;
    }
274

Dennis Gläser's avatar
Dennis Gläser committed
275
276
277
278
279
280
281
282
283
284
    //! clears all the containers
    void reset_()
    {
        elements_.clear();
        scvs_.clear();
        scvfs_.clear();
        globalLocalScvfPairedData_.clear();
        outsideScvfIndices_.clear();
        dirichletData_.clear();
    }
285

Dennis Gläser's avatar
Dennis Gläser committed
286
287
288
289
290
    //! sets up the local scvs and scvfs etc.. using the iv index set
    void setLocalScope_(DataHandle& dataHandle)
    {
        //! clear previous data
        reset_();
291

Dennis Gläser's avatar
Dennis Gläser committed
292
293
        //! number of interaction-volume-local faces
        numFaces_ = indexSet().numFaces();
294

Dennis Gläser's avatar
Dennis Gläser committed
295
        //! number of interaction-volume-local (and node-local) scvs
296
        const auto& scvIndices = indexSet().globalScvIndices();
Dennis Gläser's avatar
Dennis Gläser committed
297
        const auto numLocalScvs = scvIndices.size();
298

Dennis Gläser's avatar
Dennis Gläser committed
299
300
        //! number of global scvfs appearing in this interaction volume
        const auto numGlobalScvfs = indexSet().nodalIndexSet().numScvfs();
301
302

        //! reserve memory for local entities
Dennis Gläser's avatar
Dennis Gläser committed
303
304
305
306
        elements_.reserve(numLocalScvs);
        scvs_.reserve(numLocalScvs);
        scvfs_.reserve(numFaces_);
        dirichletData_.reserve(numFaces_);
307
308
309
        globalLocalScvfPairedData_.reserve(numGlobalScvfs);
        globalScvfIndices_.reserve(numGlobalScvfs);

Dennis Gläser's avatar
Dennis Gläser committed
310
311
312
        // store outside scvf idx set on surface grids
        if (dim < dimWorld)
            outsideScvfIndices_.reserve(numGlobalScvfs);
313

314
        // set up quantities related to sub-control volumes
Dennis Gläser's avatar
Dennis Gläser committed
315
        for (LocalIndexType scvIdxLocal = 0; scvIdxLocal < numLocalScvs; scvIdxLocal++)
316
        {
Dennis Gläser's avatar
Dennis Gläser committed
317
318
319
            const auto scvIdxGlobal = scvIndices[scvIdxLocal];
            scvs_.emplace_back(fvGeometry(), fvGeometry().scv(scvIdxGlobal), scvIdxLocal, indexSet());
            elements_.emplace_back(fvGeometry().fvGridGeometry().element(scvIdxGlobal));
320
321
        }

Dennis Gläser's avatar
Dennis Gläser committed
322
323
324
325
        // keep track of the number of unknowns etc
        numUnknowns_ = 0;
        numPotentials_ = numLocalScvs;

326
        // set up quantitites related to sub-control volume faces
Dennis Gläser's avatar
Dennis Gläser committed
327
        for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numFaces_; ++faceIdxLocal)
328
        {
Dennis Gläser's avatar
Dennis Gläser committed
329
            const auto scvfIdxGlobal = indexSet().scvfIdxGlobal(faceIdxLocal);
330
            const auto& neighborScvIndicesLocal = indexSet().neighboringLocalScvIndices(faceIdxLocal);
Dennis Gläser's avatar
Dennis Gläser committed
331
            const auto insideLocalScvIdx = neighborScvIndicesLocal[0];
332

333
            // we have to use the "inside" scv face here
Dennis Gläser's avatar
Dennis Gläser committed
334
            const auto& scvf = fvGeometry().scvf(scvfIdxGlobal);
335

Dennis Gläser's avatar
Dennis Gläser committed
336
337
            // create global/local face data for this face and the global/local map
            globalLocalScvfPairedData_.emplace_back(&scvf, LocalFaceData(faceIdxLocal, insideLocalScvIdx, false));
338
339
            globalScvfIndices_.push_back(scvf.index());

Dennis Gläser's avatar
Dennis Gläser committed
340
341
            // create iv-local scvf
            if (scvf.boundary())
342
            {
Dennis Gläser's avatar
Dennis Gläser committed
343
344
345
346
                const auto insideElement = elements_[insideLocalScvIdx];
                const auto bcTypes = problem().boundaryTypes(insideElement, scvf);

                if (bcTypes.hasOnlyDirichlet())
347
                {
Dennis Gläser's avatar
Dennis Gläser committed
348
349
                    scvfs_.emplace_back(scvf, neighborScvIndicesLocal, /*isDirichlet*/true, numPotentials_++);
                    dirichletData_.emplace_back(scvf.outsideScvIdx(), scvf.ipGlobal());
350
                }
Dennis Gläser's avatar
Dennis Gläser committed
351
352
                else
                    scvfs_.emplace_back(scvf, neighborScvIndicesLocal, /*isDirichlet*/false, numUnknowns_++);
353
            }
Dennis Gläser's avatar
Dennis Gläser committed
354
            else
355
            {
Dennis Gläser's avatar
Dennis Gläser committed
356
357
358
359
360
361
362
363
364
365
                scvfs_.emplace_back(scvf, neighborScvIndicesLocal, /*isDirichlet*/false, numUnknowns_++);

                // add outside faces to the global/local face data
                for (unsigned int i = 1; i < neighborScvIndicesLocal.size(); ++i)
                {
                    const auto outsideLocalScvIdx = neighborScvIndicesLocal[i];

                    // loop over scvfs in outside scv until we find the one coinciding with current scvf
                    for (int coord = 0; coord < dim; ++coord)
                    {
366
                        if (indexSet().scvfIdxLocal(outsideLocalScvIdx, coord) == faceIdxLocal)
Dennis Gläser's avatar
Dennis Gläser committed
367
                        {
368
                            const auto globalScvfIdx = indexSet().nodalIndexSet().scvfIdxGlobal(outsideLocalScvIdx, coord);
Dennis Gläser's avatar
Dennis Gläser committed
369
370
371
372
373
374
375
376
                            const auto& flipScvf = fvGeometry().scvf(globalScvfIdx);
                            globalLocalScvfPairedData_.emplace_back(&flipScvf, LocalFaceData(faceIdxLocal, outsideLocalScvIdx, true));
                            globalScvfIndices_.push_back(flipScvf.index());
                            if (dim < dimWorld)
                                outsideScvfIndices_.push_back(flipScvf.index());
                        }
                    }
                }
377
            }
378
        }
Dennis Gläser's avatar
Dennis Gläser committed
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393

        // resize the matrices in the data handle
        dataHandle.resizeT(numFaces_, numPotentials_);
        if (requireABMatrix_())
            dataHandle.resizeAB(numUnknowns_, numPotentials_);

        // resize the local matrices
        A_.resize(numUnknowns_, numUnknowns_);
        B_.resize(numUnknowns_, numPotentials_);
        C_.resize(numFaces_, numUnknowns_);
        D_.resize(numFaces_, numPotentials_);

        // on surface grids, resize the vector containing the "outside" transmissibilities
        if (dim < dimWorld)
            dataHandle.resizeOutsideTij(outsideScvfIndices_.size(), numPotentials_);
394
395
    }

Dennis Gläser's avatar
Dennis Gläser committed
396
    //! Assembles the local matrices that define the local system of equations and flux expressions
397
    template<typename GetTensorFunction>
Dennis Gläser's avatar
Dennis Gläser committed
398
    void assembleLocalMatrices_(const GetTensorFunction& getTensor)
399
    {
Dennis Gläser's avatar
Dennis Gläser committed
400
401
402
403
404
        // reset matrices
        A_ = 0.0;
        B_ = 0.0;
        C_ = 0.0;
        D_ = 0.0;
405

406
        // reserve space for the omegas
Dennis Gläser's avatar
Dennis Gläser committed
407
        wijk_.resize(scvfs_.size());
408

409
        // loop over the local faces
Dennis Gläser's avatar
Dennis Gläser committed
410
        for (unsigned int faceIdx = 0; faceIdx < numFaces_; ++faceIdx)
411
        {
Dennis Gläser's avatar
Dennis Gläser committed
412
            const auto& curLocalScvf = localScvf(faceIdx);
413
            const auto& curGlobalScvf = fvGeometry().scvf(curLocalScvf.globalScvfIndex());
Dennis Gläser's avatar
Dennis Gläser committed
414
415
            const auto curIsDirichlet = curLocalScvf.isDirichlet();
            const auto curLocalDofIdx = curLocalScvf.localDofIndex();
416
417

            // get diffusion tensor in "positive" sub volume
418
            const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
Dennis Gläser's avatar
Dennis Gläser committed
419
420
            const auto posLocalScvIdx = neighborScvIndices[0];
            const auto& posLocalScv = localScv(posLocalScvIdx);
421
            const auto& posGlobalScv = fvGeometry().scv(posLocalScv.globalScvIndex());
Dennis Gläser's avatar
Dennis Gläser committed
422
423
424
            const auto& posVolVars = elemVolVars()[posGlobalScv];
            const auto& posElement = element(posLocalScvIdx);
            const auto tensor = getTensor(problem(), posElement, posVolVars, fvGeometry(), posGlobalScv);
425
426

            // the omega factors of the "positive" sub volume
Dennis Gläser's avatar
Dennis Gläser committed
427
            auto posWijk = calculateOmegas_(posLocalScv, curGlobalScvf.unitOuterNormal(), curGlobalScvf.area(), tensor);
428
            posWijk *= posVolVars.extrusionFactor();
429

Dennis Gläser's avatar
Dennis Gläser committed
430
431
            // go over the coordinate directions in the positive sub volume
            for (unsigned int localDir = 0; localDir < dim; localDir++)
432
            {
433
                const auto otherLocalScvfIdx = posLocalScv.scvfIdxLocal(localDir);
Dennis Gläser's avatar
Dennis Gläser committed
434
435
                const auto& otherLocalScvf = localScvf(otherLocalScvfIdx);
                const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
436

Dennis Gläser's avatar
Dennis Gläser committed
437
438
439
                // if we are not on a Dirichlet face, add entries associated with unknown face pressures
                // i.e. in matrix C and maybe A (if current face is not a Dirichlet face)
                if (!otherLocalScvf.isDirichlet())
440
                {
Dennis Gläser's avatar
Dennis Gläser committed
441
442
443
                    C_[faceIdx][otherLocalDofIdx] -= posWijk[localDir];
                    if (!curIsDirichlet)
                        A_[curLocalDofIdx][otherLocalDofIdx] -= posWijk[localDir];
444
                }
Dennis Gläser's avatar
Dennis Gläser committed
445
446
                // the current face is a Dirichlet face and creates entries in D & maybe B
                else
447
                {
Dennis Gläser's avatar
Dennis Gläser committed
448
449
450
                    D_[faceIdx][otherLocalDofIdx] -= posWijk[localDir];
                    if (!curIsDirichlet)
                        B_[curLocalDofIdx][otherLocalDofIdx] += posWijk[localDir];
451
                }
452
453

                // add entries related to pressures at the scv centers (dofs)
Dennis Gläser's avatar
Dennis Gläser committed
454
455
                const auto posScvLocalDofIdx = posLocalScv.localDofIndex();
                D_[faceIdx][posScvLocalDofIdx] += posWijk[localDir];
456

Dennis Gläser's avatar
Dennis Gläser committed
457
458
                if (!curIsDirichlet)
                    B_[curLocalDofIdx][posScvLocalDofIdx] -= posWijk[localDir];
459
460
            }

461
            // store the omegas
Dennis Gläser's avatar
Dennis Gläser committed
462
            wijk_[faceIdx].emplace_back(std::move(posWijk));
463

Dennis Gläser's avatar
Dennis Gläser committed
464
465
            // If we are on an interior face, add values from negative sub volume
            if (!curGlobalScvf.boundary())
466
            {
467
                // loop over all the outside neighbors of this face and add entries
Dennis Gläser's avatar
Dennis Gläser committed
468
                for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside)
469
                {
Dennis Gläser's avatar
Dennis Gläser committed
470
471
                    const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1];
                    const auto& negLocalScv = localScv(negLocalScvIdx);
472
                    const auto& negGlobalScv = fvGeometry().scv(negLocalScv.globalScvIndex());
Dennis Gläser's avatar
Dennis Gläser committed
473
474
475
                    const auto& negVolVars = elemVolVars()[negGlobalScv];
                    const auto& negElement = element(negLocalScvIdx);
                    const auto negTensor = getTensor(problem(), negElement, negVolVars, fvGeometry(), negGlobalScv);
476

477
                    // the omega factors of the "negative" sub volume
478
479
480
481
482
                    DimVector negWijk;

                    // if dim < dimWorld, use outside normal vector
                    if (dim < dimWorld)
                    {
Dennis Gläser's avatar
Dennis Gläser committed
483
484
                        const auto& flipScvf = fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside);
                        auto negNormal = flipScvf.unitOuterNormal();
485
                        negNormal *= -1.0;
Dennis Gläser's avatar
Dennis Gläser committed
486
                        negWijk = calculateOmegas_(negLocalScv, negNormal, curGlobalScvf.area(), negTensor);
487
488
                    }
                    else
Dennis Gläser's avatar
Dennis Gläser committed
489
                        negWijk = calculateOmegas_(negLocalScv, curGlobalScvf.unitOuterNormal(), curGlobalScvf.area(), negTensor);
490
491

                    // scale by extrusion factpr
492
                    negWijk *= negVolVars.extrusionFactor();
493

Dennis Gläser's avatar
Dennis Gläser committed
494
                    // go over the coordinate directions in the positive sub volume
495
                    for (int localDir = 0; localDir < dim; localDir++)
496
                    {
497
                        const auto otherLocalScvfIdx = negLocalScv.scvfIdxLocal(localDir);
Dennis Gläser's avatar
Dennis Gläser committed
498
499
                        const auto& otherLocalScvf = localScvf(otherLocalScvfIdx);
                        const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
500

Dennis Gläser's avatar
Dennis Gläser committed
501
502
503
504
                        if (!otherLocalScvf.isDirichlet())
                            A_[curLocalDofIdx][otherLocalDofIdx] += negWijk[localDir];
                        else
                            B_[curLocalDofIdx][otherLocalDofIdx] -= negWijk[localDir];
505
506

                        // add entries to matrix B
Dennis Gläser's avatar
Dennis Gläser committed
507
                        B_[curLocalDofIdx][negLocalScv.localDofIndex()] += negWijk[localDir];
508
509
                    }

Dennis Gläser's avatar
Dennis Gläser committed
510
511
                    // store the omegas
                    wijk_[faceIdx].emplace_back(std::move(negWijk));
512
513
514
515
516
                }
            }
        }
    }

Dennis Gläser's avatar
Dennis Gläser committed
517
518
    //! for interaction volumes that have only dirichlet scvfs,
    //! the transmissibility matrix can be assembled directly
519
    template<typename GetTensorFunction>
Dennis Gläser's avatar
Dennis Gläser committed
520
    void assemblePureDirichletSystem_(const GetTensorFunction& getTensor, Matrix& T)
521
    {
Dennis Gläser's avatar
Dennis Gläser committed
522
523
        // reset the transmissibility matrix beforehand
        T = 0.0;
524

525
        // Loop over all the faces, in this case these are all dirichlet boundaries
Dennis Gläser's avatar
Dennis Gläser committed
526
        for (unsigned int faceIdx = 0; faceIdx < numFaces_; ++faceIdx)
527
        {
Dennis Gläser's avatar
Dennis Gläser committed
528
            const auto& curLocalScvf = localScvf(faceIdx);
529
            const auto& curGlobalScvf = fvGeometry().scvf(curLocalScvf.globalScvfIndex());
Dennis Gläser's avatar
Dennis Gläser committed
530

531
            // get diffusion tensor in "positive" sub volume
532
            const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
Dennis Gläser's avatar
Dennis Gläser committed
533
534
            const auto posLocalScvIdx = neighborScvIndices[0];
            const auto& posLocalScv = localScv(posLocalScvIdx);
535
            const auto& posGlobalScv = fvGeometry().scv(posLocalScv.globalScvIndex());
Dennis Gläser's avatar
Dennis Gläser committed
536
537
538
            const auto& posVolVars = elemVolVars()[posGlobalScv];
            const auto& posElement = element(posLocalScvIdx);
            const auto tensor = getTensor(problem(), posElement, posVolVars, fvGeometry(), posGlobalScv);
539
540

            // the omega factors of the "positive" sub volume
Dennis Gläser's avatar
Dennis Gläser committed
541
            auto posWijk = calculateOmegas_(posLocalScv, curGlobalScvf.unitOuterNormal(), curGlobalScvf.area(), tensor);
542
            posWijk *= posVolVars.extrusionFactor();
543

Dennis Gläser's avatar
Dennis Gläser committed
544
            const auto posScvLocalDofIdx = posLocalScv.localDofIndex();
545
546
            for (int localDir = 0; localDir < dim; localDir++)
            {
547
                const auto otherLocalScvfIdx = posLocalScv.scvfIdxLocal(localDir);
Dennis Gläser's avatar
Dennis Gläser committed
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
                const auto& otherLocalScvf = localScvf(otherLocalScvfIdx);
                const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
                T[faceIdx][otherLocalDofIdx] -= posWijk[localDir];
                T[faceIdx][posScvLocalDofIdx] += posWijk[localDir];
            }
        }
    }

    //! computes the transmissibilities associated with "outside" faces on surface grids
    void computeOutsideTransmissibilities_(DataHandle& dataHandle) const
    {
        assert(dim < dimWorld && "transmissibilities for outside scvfs can only be computed for dim < dimWorld!");

        for (const auto& globalLocalData : globalLocalScvfPairedData_)
        {
            const auto& localFaceData = globalLocalData.second;

            // continue only for "outside" faces
            if (!localFaceData.isOutside)
                continue;

            const auto localScvIdx = localFaceData.localScvIndex;
            const auto localScvfIdx = localFaceData.localScvfIndex;
            const auto& posLocalScv = localScv(localScvIdx);

573
            const auto idxInNeighbors = this->findIndexInVector(localScvf(localScvfIdx).neighboringLocalScvIndices(), localScvIdx);
Dennis Gläser's avatar
Dennis Gläser committed
574
575
576
577
578
579
580
581
582
            const auto& wijk = wijk_[localScvfIdx][idxInNeighbors];

            // store the calculated transmissibilities in the data handle
            const auto idxInOutsideFaces = this->findIndexInVector(outsideScvfIndices_, globalLocalData.first->index());
            auto& tij = dataHandle.outsideTij()[idxInOutsideFaces];
            tij = 0.0;

            for (int localDir = 0; localDir < dim; localDir++)
            {
583
                const auto curLocalScvfIdx = posLocalScv.scvfIdxLocal(localDir);
Dennis Gläser's avatar
Dennis Gläser committed
584
585
                const auto& curLocalScvf = localScvf(curLocalScvfIdx);
                if (!curLocalScvf.isDirichlet())
586
                {
Dennis Gläser's avatar
Dennis Gläser committed
587
588
589
                    auto tmp = dataHandle.AB()[curLocalScvf.localDofIndex()];
                    tmp *= wijk[localDir];
                    tij -= tmp;
590
591
                }
                else
Dennis Gläser's avatar
Dennis Gläser committed
592
                    tij[curLocalScvf.localDofIndex()] -= wijk[localDir];
593

Dennis Gläser's avatar
Dennis Gläser committed
594
595
                // add entry from the scv unknown
                tij[localScvIdx] += wijk[localDir];
596
597
598
599
            }
        }
    }

Dennis Gläser's avatar
Dennis Gläser committed
600
    // calculates n_i^T*K_j*nu_k
601
    DimVector calculateOmegas_(const LocalScvType& localScv,
602
603
                               const GlobalPosition normal,
                               const Scalar area,
604
605
                               const Tensor& T) const
    {
Dennis Gläser's avatar
Dennis Gläser committed
606
607
608
        // make sure we have positive definite diffsion tensors
        assert(this->tensorIsPositiveDefinite(T) && "only positive definite tensors can be handled by mpfa methods");

609
        DimVector wijk;
610
        GlobalPosition tmp;
611
612
613
        for (int dir = 0; dir < dim; ++dir)
        {
            T.mv(localScv.innerNormal(dir), tmp);
614
            wijk[dir] = tmp*normal;
615
        }
616
        wijk *= area;
617
618
619
620
621
        wijk /= localScv.detX();

        return wijk;
    }

622
    // calculates n_i^T*K_j*nu_k
623
    DimVector calculateOmegas_(const LocalScvType& localScv,
624
625
                               const GlobalPosition normal,
                               const Scalar area,
626
627
                               const Scalar t) const
    {
628
629
630
        // make sure we have positive diffusion coefficients
        assert(t > 0.0 && "non-positive diffusion coefficients cannot be handled by mpfa methods");

631
        DimVector wijk;
632
        GlobalPosition tmp(normal);
633
634
635
636
        tmp *= t;

        for (int dir = 0; dir < dim; ++dir)
            wijk[dir] = tmp*localScv.innerNormal(dir);
637
        wijk *= area;
638
639
640
641
642
        wijk /= localScv.detX();

        return wijk;
    }

643
644
645
    const Problem* problemPtr_;
    const FVElementGeometry* fvGeometryPtr_;
    const ElementVolumeVariables* elemVolVarsPtr_;
Dennis Gläser's avatar
Dennis Gläser committed
646
    const IndexSet* indexSetPtr_;
647

648
    // Variables defining the local scope
Dennis Gläser's avatar
Dennis Gläser committed
649
650
651
    std::vector<Element> elements_;
    std::vector<LocalScvType> scvs_;
    std::vector<LocalScvfType> scvfs_;
652
    std::vector<GlobalLocalFaceDataPair> globalLocalScvfPairedData_;
Dennis Gläser's avatar
Dennis Gläser committed
653
654
655
    DirichletDataContainer dirichletData_;
    GlobalIndexContainer globalScvfIndices_;
    GlobalIndexContainer outsideScvfIndices_;
656

Dennis Gläser's avatar
Dennis Gläser committed
657
658
659
660
    // sizes involved in the local matrices
    unsigned int numFaces_;
    unsigned int numUnknowns_;
    unsigned int numPotentials_;
661

Dennis Gläser's avatar
Dennis Gläser committed
662
663
    // The omega factors and the matrix A⁻1*B are stored
    // in order to recover the transmissibilities of outside faces on network grids
664
    std::vector< std::vector< DimVector > > wijk_;
Dennis Gläser's avatar
Dennis Gläser committed
665
    Matrix CAinv_;
666

Dennis Gläser's avatar
Dennis Gläser committed
667
668
669
670
671
    // Matrices involved in transmissibility calculations
    Matrix A_;
    Matrix B_;
    Matrix C_;
    Matrix D_;
672
673
674
675
676
};

} // end namespace

#endif