interactionvolume.hh 23.9 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
30
31
32
33
34
35
36
37
#include <dumux/discretization/cellcentered/mpfa/interactionvolumeseed.hh>
#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
#include <dumux/discretization/cellcentered/mpfa/methods.hh>

#include "localsubcontrolentityseeds.hh"
#include "localsubcontrolentities.hh"

namespace Dumux
{
38
//! Specialization of the interaction volume traits clas
39
40
template<class TypeTag>
class CCMpfaOInteractionVolumeTraits
41
42
43
44
{
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);

45
    static const int dim = GridView::dimension;
46
47
    static const int dimWorld = GridView::dimensionworld;
public:
48
49
    using BoundaryInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>;

50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
    using DimVector = Dune::FieldVector<Scalar, dim>;
    using Tensor = Dune::FieldMatrix<Scalar, dim, dim>;

    using LocalIndexType = std::uint8_t;
    using LocalIndexSet = std::vector<LocalIndexType>;
    using LocalIndexPair = std::pair<LocalIndexType, bool>;
    using GlobalIndexType = typename GridView::IndexSet::IndexType;
    using GlobalIndexSet = std::vector<GlobalIndexType>;

    using Matrix = Dune::DynamicMatrix<Scalar>;
    using Vector = typename Matrix::row_type;

    using PositionVector = std::vector<GlobalPosition>;

    using LocalScvSeedType = CCMpfaOLocalScvSeed<GlobalIndexType, LocalIndexType>;
    using LocalScvfSeedType = CCMpfaOLocalScvfSeed<GlobalIndexType, LocalIndexType>;
    using Seed = CCMpfaInteractionVolumeSeed<LocalScvSeedType, LocalScvfSeedType>;
};
69
70
71
72
73
/*!
 * \ingroup Mpfa
 * \brief Base class for the interaction volumes of the mpfa-o method
 */
template<class TypeTag>
74
class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> : public CCMpfaInteractionVolumeBase<TypeTag, CCMpfaOInteractionVolumeTraits<TypeTag>>
75
76
77
78
79
80
81
82
83
{
    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);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);

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

87
88
    // choose names different from those in the base class
    // base class typedefs are public to be used by other classes
89
    using Traits = CCMpfaOInteractionVolumeTraits<TypeTag>;
90
91
92
93
94
95
96
97
    using GlobalIdxType = typename Traits::GlobalIndexType;
    using LocalIdxType = typename Traits::LocalIndexType;
    using GlobalIdxSet = typename Traits::GlobalIndexSet;
    using LocalIdxSet = typename Traits::LocalIndexSet;
    using LocalIdxPair = typename Traits::LocalIndexPair;
    using IVSeed = typename Traits::Seed;

    using GlobalPosition = typename Traits::GlobalPosition;
98
99
100
    using PosVector = typename Traits::PositionVector;
    using DynamicVector = typename Traits::Vector;
    using DynamicMatrix = typename Traits::Matrix;
101
    using Tensor = typename Traits::Tensor;
102
103
104
105
106

    using LocalScvType = CCMpfaOLocalScv<TypeTag>;
    using LocalScvfType = CCMpfaOLocalScvf<TypeTag>;

public:
107
    CCMpfaInteractionVolumeImplementation(const IVSeed& seed,
108
109
110
111
112
113
                                          const Problem& problem,
                                          const FVElementGeometry& fvGeometry,
                                          const ElementVolumeVariables& elemVolVars)
    : problemRef_(problem),
      fvGeometryRef_(fvGeometry),
      elemVolVarsRef_(elemVolVars),
114
115
      onBoundary_(seed.onBoundary()),
      hasInteriorBoundary_(false)
116
117
118
119
    {
        // create local sub control entities from the seed
        createLocalEntities_(seed);

120
121
122
123
124
        // fill info on the dof stencil
        const auto numLocalScvs = localScvs_.size();
        stencil_.reserve(numLocalScvs);
        volVarsStencil_.reserve(numLocalScvs);
        volVarsPositions_.reserve(numLocalScvs);
125
126
        for (const auto& localScv : localScvs_)
        {
127
128
129
            const auto globalScvIdx = localScv.globalIndex();
            stencil_.push_back(globalScvIdx);
            volVarsStencil_.push_back(globalScvIdx);
130
131
132
            volVarsPositions_.push_back(localScv.center());
        }

133
134
        // eventually add dirichlet vol var indices and set up local index sets of flux and dirichlet faces
        LocalIdxType localScvfIdx = 0;
135
        for (const auto& localScvf : localScvfs_)
136
        {
137
138
139
140
            auto faceType = localScvf.faceType();

            // eventually add vol var index and corresponding position
            if (faceType == MpfaFaceTypes::dirichlet || faceType == MpfaFaceTypes::interiorDirichlet)
141
142
143
144
145
            {
                volVarsStencil_.push_back(localScvf.outsideGlobalScvIndex());
                volVarsPositions_.push_back(localScvf.ip());
            }

146
            // fill local scvf type index sets accordingly
147
148
149
150
            if (faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet)
                fluxFaceIndexSet_.push_back(localScvfIdx++);
            else
                dirichletFaceIndexSet_.push_back(localScvfIdx++);
151
152
153
154

            // if face is on an interior boundary, save this information
            if (!hasInteriorBoundary_ && (faceType == MpfaFaceTypes::interiorDirichlet || faceType == MpfaFaceTypes::interiorNeumann))
                hasInteriorBoundary_ = true;
155
156
        }

157
        neumannFluxes_ = DynamicVector(fluxFaceIndexSet_.size(), 0.0);
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
    }

    template<typename GetTensorFunction>
    void solveLocalSystem(const GetTensorFunction& getTensor)
    {
        const std::size_t numFluxFaces = fluxScvfIndexSet_().size();

        // if only dirichlet faces are present, assemble T_ directly
        if (numFluxFaces == 0)
            return assemblePureDirichletSystem_(getTensor);

        const std::size_t numFaces = localScvfs_.size();
        const std::size_t numPotentials = volVarsStencil().size();

        // the local matrices
173
174
175
176
        DynamicMatrix A(numFluxFaces, numFluxFaces, 0.0);
        DynamicMatrix B(numFluxFaces, numPotentials, 0.0);
        DynamicMatrix C(numFaces, numFluxFaces, 0.0);
        DynamicMatrix D(numFaces, numPotentials, 0.0);
177
178
179
180

        assembleLocalMatrices_(getTensor, A, B, C, D);

        // solve local system and store matrices
181
        DynamicMatrix copy(B);
182
183
184
185
186
187
188
189
        A.invert();
        AinvB_ = B.leftmultiply(A);
        CAinv_ = C.rightmultiply(A);
        T_ = multiplyMatrices(CAinv_, copy);
        T_ += D;
    }

    template<typename UpwindFactorFunction>
190
    void assembleNeumannFluxes(const UpwindFactorFunction& upwindFactor, const unsigned int eqIdx)
191
    {
192
193
        if (!onBoundary())
            return;
194

195
        LocalIdxType fluxFaceIdx = 0;
196
197
198
199
200
201
        for (const auto localFluxFaceIdx : fluxFaceIndexSet_)
        {
            const auto& localScvf = localScvf_(localFluxFaceIdx);
            const auto faceType = localScvf.faceType();
            if (faceType == MpfaFaceTypes::neumann || faceType == MpfaFaceTypes::interiorNeumann)
            {
202
203
204
205
                // boundary neumann fluxes stay zero when tpfa boundary handling is on
                if (GET_PROP_VALUE(TypeTag, UseTpfaBoundary) && faceType == MpfaFaceTypes::neumann)
                    continue;

206
207
208
                const auto& element = localElement_(localScvf.insideLocalScvIndex());
                const auto& globalScvf = fvGeometry_().scvf(localScvf.insideGlobalScvfIndex());
                auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf)[eqIdx];
209
                neumannFlux *= globalScvf.area();
210

211
212
213
214
215
216
217
                // if we are not on an interior neumann boundary, we have to recover -k*gradh
                if (faceType == MpfaFaceTypes::neumann)
                {
                    const auto& insideScv = fvGeometry_().scv(globalScvf.insideScvIdx());
                    const auto& volVars = elemVolVars_()[insideScv];
                    neumannFlux /= upwindFactor(volVars);
                }
218
219
220
221
222
223
224
225

                neumannFluxes_[fluxFaceIdx] = neumannFlux;
            }

            fluxFaceIdx++;
        }
    }

226
    bool onBoundary() const
227
228
    { return onBoundary_; }

229
230
231
232
    bool hasInteriorBoundary() const
    { return hasInteriorBoundary_; }

    const GlobalIdxSet& stencil() const
233
234
    { return stencil_; }

235
    const GlobalIdxSet& volVarsStencil() const
236
237
    { return volVarsStencil_; }

238
    const PosVector& volVarsPositions() const
239
240
    { return volVarsPositions_; }

241
    GlobalIdxSet globalScvfs() const
242
    {
243
244
        GlobalIdxSet globalScvfs;
        globalScvfs.reserve(localScvfs_.size());
245
246
247

        for (const auto& localScvf : localScvfs_)
        {
248
            globalScvfs.push_back(localScvf.insideGlobalScvfIndex());
249
            if (!localScvf.boundary())
250
251
252
253
254
255
256
257
                globalScvfs.push_back(localScvf.outsideGlobalScvfIndex());
        }

        // if interior boundaries are present we have to make the entries unique
        if (hasInteriorBoundary())
        {
            std::sort(globalScvfs.begin(), globalScvfs.end());
            globalScvfs.erase(std::unique(globalScvfs.begin(), globalScvfs.end()), globalScvfs.end());
258
259
        }

260
        return globalScvfs;
261
262
    }

263
    LocalIdxPair getLocalIndexPair(const SubControlVolumeFace& scvf) const
264
265
266
    {
        auto scvfGlobalIdx = scvf.index();

267
        LocalIdxType localIdx = 0;
268
269
270
271
272
273
274
275
276
277
278
279
280
281
        for (const auto& localScvf : localScvfs_)
        {
            if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
                return std::make_pair(localIdx, false);

            if (!localScvf.boundary() && localScvf.outsideGlobalScvfIndex() == scvfGlobalIdx)
                return std::make_pair(localIdx, true);

            localIdx++;
        }

        DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv face in the interaction volume for the given scvf with index: " << scvf.index());
    }

282
    DynamicVector getTransmissibilities(const std::pair<LocalIdxType, bool>& localIndexPair) const
283
284
285
286
287
288
289
290
    {
        auto tij = T_[localIndexPair.first];

        if (localIndexPair.second)
            tij *= -1.0;
        return tij;
    }

291
    Scalar getNeumannFlux(const std::pair<LocalIdxType, bool>& localIndexPair) const
292
    {
293
294
        if (fluxScvfIndexSet_().size() == 0)
            return 0.0;
295

296
        auto flux = CAinv_[localIndexPair.first] * neumannFluxes_;
297
298
299
300
301
302
303
        if (localIndexPair.second)
            return -1.0*flux;
        return flux;
    }

private:

304
    const LocalScvfType& localScvf_(const LocalIdxType localScvfIdx) const
305
306
    { return localScvfs_[localScvfIdx]; }

307
    const LocalScvType& localScv_(const LocalIdxType localScvIdx) const
308
309
    { return localScvs_[localScvIdx]; }

310
    const LocalIdxSet& fluxScvfIndexSet_() const
311
312
    { return fluxFaceIndexSet_; }

313
    const LocalIdxSet& dirichletScvfIndexSet_() const
314
315
    { return dirichletFaceIndexSet_; }

316
    const Element& localElement_(const LocalIdxType localScvIdx) const
317
318
    { return localElements_[localScvIdx]; }

319
    void createLocalEntities_(const IVSeed& seed)
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
    {
        auto numScvs = seed.scvSeeds().size();
        auto numScvfs = seed.scvfSeeds().size();

        localElements_.reserve(numScvs);
        localScvs_.reserve(numScvs);
        localScvfs_.reserve(numScvfs);

        for (const auto& scvSeed : seed.scvSeeds())
        {
            auto element = problem_().model().globalFvGeometry().element(scvSeed.globalIndex());
            localScvs_.emplace_back(LocalScvType(element, fvGeometry_(), scvSeed));
            localElements_.emplace_back(std::move(element));
        }

        for (const auto& scvfSeed : seed.scvfSeeds())
        {
            // we have to use the "inside" scv face here
            const auto& scvf = fvGeometry_().scvf(scvfSeed.insideGlobalScvfIndex());
            localScvfs_.emplace_back(LocalScvfType(scvfSeed, scvf));
        }
    }

    template<typename GetTensorFunction>
    void assembleLocalMatrices_(const GetTensorFunction& getTensor,
345
346
347
348
                                DynamicMatrix& A,
                                DynamicMatrix& B,
                                DynamicMatrix& C,
                                DynamicMatrix& D)
349
350
351
352
353
    {
        const Scalar xi = GET_PROP_VALUE(TypeTag, Xi);
        const std::size_t numLocalScvs = localScvs_.size();

        // loop over the local faces
354
        LocalIdxType rowIdx = 0;
355
356
357
358
        for (const auto& localScvf : localScvfs_)
        {
            auto faceType = localScvf.faceType();
            bool hasUnknown = faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet;
359
            LocalIdxType idxInFluxFaces = hasUnknown ? this->findLocalIndex(fluxScvfIndexSet_(), rowIdx) : -1;
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383

            // get diffusion tensor in "positive" sub volume
            const auto posLocalScvIdx = localScvf.insideLocalScvIndex();
            const auto& posLocalScv = localScv_(posLocalScvIdx);
            const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
            auto element = localElement_(posLocalScvIdx);
            auto tensor = getTensor(element, elemVolVars_()[posGlobalScv], posGlobalScv);

            // the omega factors of the "positive" sub volume
            auto posWijk = calculateOmegas_(posLocalScv, localScvf, tensor);
            posWijk *= problem_().boxExtrusionFactor(element, posGlobalScv);

            // Check the local directions of the positive sub volume
            for (int localDir = 0; localDir < dim; localDir++)
            {
                auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir);
                const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
                auto curFaceType = curLocalScvf.faceType();
                bool curFaceHasUnknown = curFaceType != MpfaFaceTypes::dirichlet && curFaceType != MpfaFaceTypes::interiorDirichlet;

                // First, add the entries associated with face pressures (unkown or dirichlet)
                if (curFaceHasUnknown)
                {
                    // we need the index of the current local scvf in the flux face indices
384
                    auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx);
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402

                    C[rowIdx][curIdxInFluxFaces] += posWijk[localDir];
                    if (hasUnknown)
                    {
                        if (faceType == MpfaFaceTypes::interiorNeumann)
                        {
                            A[idxInFluxFaces][curIdxInFluxFaces] += xi*posWijk[localDir];
                            // If facet coupling active, eventually add entry coming from the lowdim domain
                            if (GET_PROP_VALUE(TypeTag, FacetCoupling) && curIdxInFluxFaces == idxInFluxFaces)
                                DUNE_THROW(Dune::NotImplemented, "Facet coupling");
                        }
                        else
                            A[idxInFluxFaces][curIdxInFluxFaces] += posWijk[localDir];
                    }
                }
                else
                {
                    // the current face is a Dirichlet face and creates entries in D & eventually B
403
                    auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
                    D[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];

                    if (hasUnknown)
                        B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] -= posWijk[localDir];
                }

                // add entries related to pressures at the scv centers (dofs)
                D[rowIdx][posLocalScvIdx] -= posWijk[localDir];

                if (hasUnknown)
                {
                    if (faceType == MpfaFaceTypes::interiorNeumann)
                        B[idxInFluxFaces][posLocalScvIdx] += xi*posWijk[localDir];
                    else
                        B[idxInFluxFaces][posLocalScvIdx] += posWijk[localDir];
                }
            }

            // If not on a boundary or interior dirichlet face, add entries for the "negative" scv
            if (faceType == MpfaFaceTypes::interior || faceType == MpfaFaceTypes::interiorNeumann)
            {
                const auto negLocalScvIdx = localScvf.outsideLocalScvIndex();
                const auto& negLocalScv = localScv_(negLocalScvIdx);
                const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
                auto negElement = localElement_(negLocalScvIdx);;
                auto negTensor = getTensor(negElement, elemVolVars_()[negGlobalScv], negGlobalScv);

                // the omega factors of the "negative" sub volume
                auto negWijk = calculateOmegas_(negLocalScv, localScvf, negTensor);
                negWijk *= problem_().boxExtrusionFactor(negElement, negGlobalScv);

                // Check local directions of negative sub volume
                for (int localDir = 0; localDir < dim; localDir++)
                {
                    auto curLocalScvfIdx = negLocalScv.localScvfIndex(localDir);
                    const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
                    auto curFaceType = curLocalScvf.faceType();
                    bool curFaceHasUnknown = curFaceType != MpfaFaceTypes::dirichlet && curFaceType != MpfaFaceTypes::interiorDirichlet;

                    if (curFaceHasUnknown)
                    {
                        // we need the index of the current local scvf in the flux face indices
446
                        auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx);
447
448
449
450
451
452
453
454
455

                        if (faceType == MpfaFaceTypes::interiorNeumann)
                            A[idxInFluxFaces][curIdxInFluxFaces] -= (1.0-xi)*negWijk[localDir];
                        else
                            A[idxInFluxFaces][curIdxInFluxFaces] -= negWijk[localDir];
                    }
                    else
                    {
                        // the current face is a Dirichlet face and creates entries in B
456
                        auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
                        B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] += negWijk[localDir];
                    }

                    // add entries to matrix B
                    if (faceType == MpfaFaceTypes::interiorNeumann)
                        B[idxInFluxFaces][negLocalScvIdx] -= (1.0-xi)*negWijk[localDir];
                    else
                        B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir];
                }
            }
            // go to the next face
            rowIdx++;
        }
    }

    template<typename GetTensorFunction>
    void assemblePureDirichletSystem_(const GetTensorFunction& getTensor)
    {
        const std::size_t numLocalScvs = localScvs_.size();
        const std::size_t numFaces = localScvfs_.size();
        const std::size_t numPotentials = volVarsStencil().size();

        // resize matrices, only T_ will have entries
        T_.resize(numFaces, numPotentials, 0.0);
        AinvB_.resize(0, 0);
        CAinv_.resize(0, 0);

        // Loop over all the faces, in this case these are all dirichlet boundaries
485
        LocalIdxType rowIdx = 0;
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
        for (const auto& localScvf : localScvfs_)
        {
            // get diffusion tensor in "positive" sub volume
            const auto posLocalScvIdx = localScvf.insideLocalScvIndex();
            const auto& posLocalScv = localScv_(posLocalScvIdx);
            const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
            auto element = localElement_(posLocalScvIdx);
            auto tensor = getTensor(element, elemVolVars_()[posGlobalScv], posGlobalScv);

            // the omega factors of the "positive" sub volume
            auto posWijk = calculateOmegas_(posLocalScv, localScvf, tensor);
            posWijk *= problem_().boxExtrusionFactor(element, posGlobalScv);

            for (int localDir = 0; localDir < dim; localDir++)
            {
                auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir);
502
                auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
503
504
505
506
507
508
509
510
511
512

                T_[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
                T_[rowIdx][posLocalScvIdx] -= posWijk[localDir];
            }

            // go to the next face
            rowIdx++;
        }
    }

513
    GlobalPosition calculateOmegas_(const LocalScvType& localScv,
514
515
516
                               const LocalScvfType& localScvf,
                               const Tensor& T) const
    {
517
518
        GlobalPosition wijk;
        GlobalPosition tmp;
519
520
521
522
523
524
525
        for (int dir = 0; dir < dim; ++dir)
        {
            T.mv(localScv.innerNormal(dir), tmp);
            wijk[dir] = tmp*localScvf.unitOuterNormal();
        }
        wijk *= localScvf.area();
        wijk /= localScv.detX();
526
        wijk *= -1.0;
527
528
529
530

        return wijk;
    }

531
    GlobalPosition calculateOmegas_(const LocalScvType& localScv,
532
533
534
                               const LocalScvfType& localScvf,
                               const Scalar t) const
    {
535
536
        GlobalPosition wijk;
        GlobalPosition tmp(localScvf.unitOuterNormal());
537
538
539
540
541
542
        tmp *= t;

        for (int dir = 0; dir < dim; ++dir)
            wijk[dir] = tmp*localScv.innerNormal(dir);
        wijk *= localScvf.area();
        wijk /= localScv.detX();
543
        wijk *= -1.0;
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561

        return wijk;
    }

    const Problem& problem_() const
    { return problemRef_; }

    const FVElementGeometry& fvGeometry_() const
    { return fvGeometryRef_; }

    const ElementVolumeVariables& elemVolVars_() const
    { return elemVolVarsRef_; }

    const Problem& problemRef_;
    const FVElementGeometry& fvGeometryRef_;
    const ElementVolumeVariables& elemVolVarsRef_;

    bool onBoundary_;
562
563
    bool hasInteriorBoundary_;
    LocalIdxType eqIdx_;
564
565
566
567
568

    std::vector<Element> localElements_;
    std::vector<LocalScvType> localScvs_;
    std::vector<LocalScvfType> localScvfs_;

569
570
571
    GlobalIdxSet globalScvfIndices_;
    GlobalIdxSet stencil_;
    GlobalIdxSet volVarsStencil_;
572
    PosVector volVarsPositions_;
573

574
575
    LocalIdxSet fluxFaceIndexSet_;
    LocalIdxSet dirichletFaceIndexSet_;
576

577
578
579
    DynamicMatrix T_;
    DynamicMatrix AinvB_;
    DynamicMatrix CAinv_;
580

581
    DynamicVector neumannFluxes_;
582
583
584
585
586
};

} // end namespace

#endif