interactionvolume.hh 21.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
#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
#include <dumux/discretization/cellcentered/mpfa/methods.hh>

32
#include "interactionvolumeseed.hh"
33
34
35
36
#include "localsubcontrolentities.hh"

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

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

49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
    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>;

64
65
66
    using LocalScvType = CCMpfaOLocalScv<TypeTag>;
    using LocalScvfType = CCMpfaOLocalScvf<TypeTag>;

67
    using Seed = CCMpfaOInteractionVolumeSeed<GlobalIndexType, LocalIndexType>;
68
};
69
70
71
72

//! Forward declaration of the mpfa-o interaction volume
template<class TypeTag, class Traits> class CCMpfaOInteractionVolume;

73
74
/*!
 * \ingroup Mpfa
75
76
77
78
 * \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.
79
80
 */
template<class TypeTag>
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> : public CCMpfaOInteractionVolume<TypeTag, CCMpfaOInteractionVolumeTraits<TypeTag>>
{
    using Traits = CCMpfaOInteractionVolumeTraits<TypeTag>;
    using ParentType = CCMpfaOInteractionVolume<TypeTag, Traits>;

    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);

    using IVSeed = typename Traits::Seed;
public:
    CCMpfaInteractionVolumeImplementation(const IVSeed& seed,
                                          const Problem& problem,
                                          const FVElementGeometry& fvGeometry,
                                          const ElementVolumeVariables& elemVolVars)
    : ParentType(seed, problem, fvGeometry, elemVolVars)
    {}

};

template<class TypeTag, class Traits>
class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase<TypeTag, Traits>
103
{
104
105
    // The interaction volume implementation has to be friend
    friend typename GET_PROP_TYPE(TypeTag, InteractionVolume);
106
    using ParentType = CCMpfaInteractionVolumeBase<TypeTag, Traits>;
107

108
109
110
111
112
113
114
115
    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);

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

119
    using DimVector = typename Traits::DimVector;
120
    using GlobalPosition = typename Traits::GlobalPosition;
121
122
    using DynamicVector = typename Traits::Vector;
    using DynamicMatrix = typename Traits::Matrix;
123
    using Tensor = typename Traits::Tensor;
124

125
126
    using LocalScvType = typename Traits::LocalScvType;
    using LocalScvfType = typename Traits::LocalScvfType;
127
128

public:
129
130
131
132
133
134
135
136
    using typename ParentType::LocalIndexType;
    using typename ParentType::LocalIndexSet;
    using typename ParentType::LocalIndexPair;
    using typename ParentType::GlobalIndexSet;
    using typename ParentType::PositionVector;
    using typename ParentType::Seed;

    CCMpfaOInteractionVolume(const Seed& seed,
137
138
139
                             const Problem& problem,
                             const FVElementGeometry& fvGeometry,
                             const ElementVolumeVariables& elemVolVars)
140
141
142
    : problemPtr_(&problem),
      fvGeometryPtr_(&fvGeometry),
      elemVolVarsPtr_(&elemVolVars),
143
144
      onBoundary_(seed.onBoundary()),
      globalScvfIndices_(seed.globalScvfIndices())
145
146
147
148
    {
        // create local sub control entities from the seed
        createLocalEntities_(seed);

149
        // reserve memory
150
        const auto numLocalScvs = localScvs_.size();
151
152
        const auto numLocalScvfs = localScvfs_.size();
        const auto maxNumVolVars = numLocalScvs + numLocalScvfs;
153
        volVarsStencil_ = seed.globalScvIndices(); // boundary vol vars are placed at the end
154
155
156
157
158
        volVarsStencil_.reserve(maxNumVolVars);
        volVarsPositions_.reserve(maxNumVolVars);
        dirichletFaceIndexSet_.reserve(numLocalScvfs);
        fluxFaceIndexSet_.reserve(numLocalScvfs);

159
        // the positions where the vol vars are defined at (required for the gravitational acceleration)
160
161
162
        for (const auto& localScv : localScvs_)
            volVarsPositions_.push_back(localScv.center());

163
        // eventually add dirichlet vol var indices and set up local index sets of flux and dirichlet faces
164
        LocalIndexType localScvfIdx = 0;
165
        for (const auto& localScvf : localScvfs_)
166
        {
167
            // eventually add vol var index and corresponding position
168
            if (localScvf.faceType() == MpfaFaceTypes::dirichlet)
169
170
171
            {
                volVarsStencil_.push_back(localScvf.outsideGlobalScvIndex());
                volVarsPositions_.push_back(localScvf.ip());
172
                dirichletFaceIndexSet_.push_back(localScvfIdx++);
173
174
            }
            else
175
                fluxFaceIndexSet_.push_back(localScvfIdx++);
176
177
        }

178
        // initialize the neumann fluxes vector to zero
179
        neumannFluxes_ = DynamicVector(fluxFaceIndexSet_.size(), 0.0);
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
    }

    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
195
196
197
198
        DynamicMatrix A(numFluxFaces, numFluxFaces, 0.0);
        DynamicMatrix B(numFluxFaces, numPotentials, 0.0);
        DynamicMatrix C(numFaces, numFluxFaces, 0.0);
        DynamicMatrix D(numFaces, numPotentials, 0.0);
199
200
201
202

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

        // solve local system and store matrices
203
        DynamicMatrix copy(B);
204
205
206
207
208
209
210
211
        A.invert();
        AinvB_ = B.leftmultiply(A);
        CAinv_ = C.rightmultiply(A);
        T_ = multiplyMatrices(CAinv_, copy);
        T_ += D;
    }

    template<typename UpwindFactorFunction>
212
    void assembleNeumannFluxes(const UpwindFactorFunction& upwindFactor, const unsigned int eqIdx)
213
    {
214
        if (!onBoundary() || GET_PROP_VALUE(TypeTag, UseTpfaBoundary))
215
            return;
216

217
        LocalIndexType fluxFaceIdx = 0;
218
219
220
221
        for (const auto localFluxFaceIdx : fluxFaceIndexSet_)
        {
            const auto& localScvf = localScvf_(localFluxFaceIdx);
            const auto faceType = localScvf.faceType();
222
            if (faceType == MpfaFaceTypes::neumann)
223
224
225
226
            {
                const auto& element = localElement_(localScvf.insideLocalScvIndex());
                const auto& globalScvf = fvGeometry_().scvf(localScvf.insideGlobalScvfIndex());
                auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf)[eqIdx];
227
                neumannFlux *= globalScvf.area();
228

229
230
231
232
                // recover -k*gradh
                const auto& insideScv = fvGeometry_().scv(globalScvf.insideScvIdx());
                const auto& volVars = elemVolVars_()[insideScv];
                neumannFlux /= upwindFactor(volVars);
233
234
235
236
237
238
239
240

                neumannFluxes_[fluxFaceIdx] = neumannFlux;
            }

            fluxFaceIdx++;
        }
    }

241
    LocalIndexPair getLocalIndexPair(const SubControlVolumeFace& scvf) const
242
243
244
    {
        auto scvfGlobalIdx = scvf.index();

245
        LocalIndexType localIdx = 0;
246
247
248
249
250
251
252
253
254
255
256
257
258
259
        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());
    }

260
    DynamicVector getTransmissibilities(const LocalIndexPair& localIndexPair) const
261
262
263
264
265
266
267
268
    {
        auto tij = T_[localIndexPair.first];

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

269
    Scalar getNeumannFlux(const LocalIndexPair& localIndexPair) const
270
    {
271
        if (!onBoundary() || fluxScvfIndexSet_().size() == 0 || GET_PROP_VALUE(TypeTag, UseTpfaBoundary))
272
            return 0.0;
273

274
        auto flux = CAinv_[localIndexPair.first] * neumannFluxes_;
275
276
277
278
279
        if (localIndexPair.second)
            return -1.0*flux;
        return flux;
    }

280
281
282
283
284
285
286
287
288
    bool onBoundary() const
    { return onBoundary_; }

    const GlobalIndexSet& volVarsStencil() const
    { return volVarsStencil_; }

    const PositionVector& volVarsPositions() const
    { return volVarsPositions_; }

289
290
291
    const GlobalIndexSet& globalScvfs() const
    { return globalScvfIndices_; }

292
293
private:

294
    const LocalScvfType& localScvf_(const LocalIndexType localScvfIdx) const
295
296
    { return localScvfs_[localScvfIdx]; }

297
    const LocalScvType& localScv_(const LocalIndexType localScvIdx) const
298
299
    { return localScvs_[localScvIdx]; }

300
    const LocalIndexSet& fluxScvfIndexSet_() const
301
302
    { return fluxFaceIndexSet_; }

303
    const LocalIndexSet& dirichletScvfIndexSet_() const
304
305
    { return dirichletFaceIndexSet_; }

306
    const Element& localElement_(const LocalIndexType localScvIdx) const
307
308
    { return localElements_[localScvIdx]; }

309
    void createLocalEntities_(const Seed& seed)
310
311
312
313
314
315
316
317
318
319
320
    {
        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());
321
            localScvs_.emplace_back(LocalScvType(problem_(), element, fvGeometry_(), scvSeed));
322
323
324
325
326
327
328
329
330
331
332
333
334
            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,
335
336
337
338
                                DynamicMatrix& A,
                                DynamicMatrix& B,
                                DynamicMatrix& C,
                                DynamicMatrix& D)
339
340
341
342
    {
        const std::size_t numLocalScvs = localScvs_.size();

        // loop over the local faces
343
        LocalIndexType rowIdx = 0;
344
345
346
        for (const auto& localScvf : localScvfs_)
        {
            auto faceType = localScvf.faceType();
347
            bool hasUnknown = faceType != MpfaFaceTypes::dirichlet;
348
            LocalIndexType idxInFluxFaces = hasUnknown ? this->findLocalIndex(fluxScvfIndexSet_(), rowIdx) : -1;
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367

            // 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);

                // First, add the entries associated with face pressures (unkown or dirichlet)
368
                if (curLocalScvf.faceType() != MpfaFaceTypes::dirichlet)
369
370
                {
                    // we need the index of the current local scvf in the flux face indices
371
                    auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx);
372
373
374

                    C[rowIdx][curIdxInFluxFaces] += posWijk[localDir];
                    if (hasUnknown)
375
                        A[idxInFluxFaces][curIdxInFluxFaces] += posWijk[localDir];
376
377
378
379
                }
                else
                {
                    // the current face is a Dirichlet face and creates entries in D & eventually B
380
                    auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
381

382
                    D[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
383
384
385
386
387
388
389
390
                    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)
391
                    B[idxInFluxFaces][posLocalScvIdx] += posWijk[localDir];
392
393
            }

394
            // If face is not on a boundary or an interior dirichlet face, add entries for the "negative" scv
395
            if (faceType == MpfaFaceTypes::interior)
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
            {
                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);

413
                    if (curLocalScvf.faceType() != MpfaFaceTypes::dirichlet)
414
415
                    {
                        // we need the index of the current local scvf in the flux face indices
416
                        auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx);
417
                        A[idxInFluxFaces][curIdxInFluxFaces] -= negWijk[localDir];
418
419
420
421
                    }
                    else
                    {
                        // the current face is a Dirichlet face and creates entries in B
422
                        auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
423
424
425
426
                        B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] += negWijk[localDir];
                    }

                    // add entries to matrix B
427
                    B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir];
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
                }
            }
            // 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
448
        LocalIndexType rowIdx = 0;
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
        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);
465
                auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
466
467
468
469
470
471
472
473
474
475

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

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

476
    DimVector calculateOmegas_(const LocalScvType& localScv,
477
478
479
                               const LocalScvfType& localScvf,
                               const Tensor& T) const
    {
480
481
        GlobalPosition wijk;
        GlobalPosition tmp;
482
483
484
485
486
487
488
        for (int dir = 0; dir < dim; ++dir)
        {
            T.mv(localScv.innerNormal(dir), tmp);
            wijk[dir] = tmp*localScvf.unitOuterNormal();
        }
        wijk *= localScvf.area();
        wijk /= localScv.detX();
489
        wijk *= -1.0;
490
491
492
493

        return wijk;
    }

494
    DimVector calculateOmegas_(const LocalScvType& localScv,
495
496
497
                               const LocalScvfType& localScvf,
                               const Scalar t) const
    {
498
499
        GlobalPosition wijk;
        GlobalPosition tmp(localScvf.unitOuterNormal());
500
501
502
503
504
505
        tmp *= t;

        for (int dir = 0; dir < dim; ++dir)
            wijk[dir] = tmp*localScv.innerNormal(dir);
        wijk *= localScvf.area();
        wijk /= localScv.detX();
506
        wijk *= -1.0;
507
508
509
510
511

        return wijk;
    }

    const Problem& problem_() const
512
    { return *problemPtr_; }
513
514

    const FVElementGeometry& fvGeometry_() const
515
    { return *fvGeometryPtr_; }
516
517

    const ElementVolumeVariables& elemVolVars_() const
518
    { return *elemVolVarsPtr_; }
519

520
521
522
    const Problem* problemPtr_;
    const FVElementGeometry* fvGeometryPtr_;
    const ElementVolumeVariables* elemVolVarsPtr_;
523
524

    bool onBoundary_;
525
    LocalIndexType eqIdx_;
526
527
528
529
530

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

531
532
533
    GlobalIndexSet globalScvfIndices_;
    GlobalIndexSet volVarsStencil_;
    PositionVector volVarsPositions_;
534

535
536
    LocalIndexSet fluxFaceIndexSet_;
    LocalIndexSet dirichletFaceIndexSet_;
537

538
539
540
    DynamicMatrix T_;
    DynamicMatrix AinvB_;
    DynamicMatrix CAinv_;
541

542
    DynamicVector neumannFluxes_;
543
544
545
546
547
};

} // end namespace

#endif