interactionvolume.hh 40.1 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
template<class TypeTag>
39
class CCMpfaOInteractionVolumeTraits : public CCMpfaInteractionVolumeTraitsBase<TypeTag>
40
{
41
    using BaseTraits = CCMpfaInteractionVolumeTraitsBase<TypeTag>;
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
48
    static const int dimWorld = GridView::dimensionworld;
    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;

49
50
public:
    using BoundaryInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>;
51

52
    using PositionVector = std::vector<GlobalPosition>;
53
54
55
    using Matrix = Dune::DynamicMatrix<Scalar>;
    using Vector = typename Matrix::row_type;

56
57
    using LocalScvType = CCMpfaOLocalScv<TypeTag>;
    using LocalScvfType = CCMpfaOLocalScvf<TypeTag>;
58
59
60
    using typename BaseTraits::LocalIndexSet;
    using typename BaseTraits::GlobalIndexSet;
    using Seed = CCMpfaOInteractionVolumeSeed<GlobalIndexSet, LocalIndexSet, dim, dimWorld>;
61
};
62
63
64
65

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

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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
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>
96
{
97
98
    // The interaction volume implementation has to be friend,
    // because some methods use the mpfa-o interaction volume as base
99
    friend typename GET_PROP_TYPE(TypeTag, InteractionVolume);
100
    using Implementation = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
101
    using ParentType = CCMpfaInteractionVolumeBase<TypeTag, Traits>;
102

103
104
105
106
    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);
107
    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
108
109
110
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
111
112
113
114
115
    using InteriorBoundaryData = typename GET_PROP_TYPE(TypeTag, InteriorBoundaryData);

    static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
    static constexpr bool facetCoupling = GET_PROP_VALUE(TypeTag, MpfaFacetCoupling);
    static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
116

117
    static const int dim = GridView::dimension;
118
    static const int dimWorld = GridView::dimensionworld;
119
    using Element = typename GridView::template Codim<0>::Entity;
120
121
    using DimVector = Dune::FieldVector<Scalar, dim>;
    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
122

123
124
    using DynamicVector = typename Traits::Vector;
    using DynamicMatrix = typename Traits::Matrix;
125
    using Tensor = typename Traits::Tensor;
126

127
128
    using LocalScvType = typename Traits::LocalScvType;
    using LocalScvfType = typename Traits::LocalScvfType;
129
130

public:
131
132
    using typename ParentType::LocalIndexType;
    using typename ParentType::LocalIndexSet;
133
    using typename ParentType::LocalFaceData;
134
135
136
137
138
    using typename ParentType::GlobalIndexSet;
    using typename ParentType::PositionVector;
    using typename ParentType::Seed;

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

151
        // reserve memory
Dennis Gläser's avatar
Dennis Gläser committed
152
153
154
        auto numLocalScvs = localScvs_.size();
        auto numLocalScvfs = localScvfs_.size();
        auto maxNumVolVars = numLocalScvs + numLocalScvfs;
155
        volVarsStencil_ = seed.globalScvIndices(); // boundary vol vars are placed at the end
156
157
158
        volVarsStencil_.reserve(maxNumVolVars);
        volVarsPositions_.reserve(maxNumVolVars);
        dirichletFaceIndexSet_.reserve(numLocalScvfs);
159
160
161
        interiorDirichletFaceIndexSet_.reserve(numLocalScvfs);
        interiorBoundaryFaceIndexSet_.reserve(numLocalScvfs);
        interiorBoundaryData_.reserve(numLocalScvfs);
162
163
        fluxFaceIndexSet_.reserve(numLocalScvfs);

164
        // the positions where the vol vars are defined at (required for the gravitational acceleration)
Dennis Gläser's avatar
Dennis Gläser committed
165
        for (auto&& localScv : localScvs_)
166
167
            volVarsPositions_.push_back(localScv.center());

168
        // maybe add dirichlet vol var indices and set up local index sets of flux and dirichlet faces
169
        LocalIndexType localScvfIdx = 0;
Dennis Gläser's avatar
Dennis Gläser committed
170
        for (auto&& localScvf : localScvfs_)
171
        {
172
173
            auto faceType = localScvf.faceType();

174
            // eventually add vol var index and corresponding position
175
            if (faceType == MpfaFaceTypes::dirichlet)
176
177
178
            {
                volVarsStencil_.push_back(localScvf.outsideGlobalScvIndex());
                volVarsPositions_.push_back(localScvf.ip());
179
                dirichletFaceIndexSet_.push_back(localScvfIdx++);
180
            }
181
182
            else if (faceType == MpfaFaceTypes::interior || faceType == MpfaFaceTypes::neumann)
            {
183
                fluxFaceIndexSet_.push_back(localScvfIdx++);
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
            }
            else if (faceType == MpfaFaceTypes::interiorNeumann)
            {
                interiorBoundaryData_.push_back(InteriorBoundaryData(problem,
                                                                     localScvf.insideGlobalScvIndex(),
                                                                     localScvf.insideGlobalScvfIndex(),
                                                                     fluxFaceIndexSet_.size(),
                                                                     faceType));
                fluxFaceIndexSet_.push_back(localScvfIdx);
                interiorBoundaryFaceIndexSet_.push_back(localScvfIdx++);
            }
            else if(faceType == MpfaFaceTypes::interiorDirichlet)
            {
                interiorBoundaryData_.push_back(InteriorBoundaryData(problem,
                                                                     localScvf.insideGlobalScvIndex(),
                                                                     localScvf.insideGlobalScvfIndex(),
                                                                     interiorDirichletFaceIndexSet_.size(),
                                                                     faceType));
                interiorDirichletFaceIndexSet_.push_back(localScvfIdx);
                interiorBoundaryFaceIndexSet_.push_back(localScvfIdx++);
            }
            else
                DUNE_THROW(Dune::InvalidStateException, "Unknown mpfa face type.");
207
208
        }

209
210
        // initialize the vector containing the neumann fluxes
        assembleNeumannFluxVector_();
211
212
213
214
215
    }

    template<typename GetTensorFunction>
    void solveLocalSystem(const GetTensorFunction& getTensor)
    {
216
        const auto numFluxFaces = fluxScvfIndexSet_().size();
217
218
219
220
221

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

222
223
        const auto numFaces = localScvfs_.size();
        const auto numPotentials = volVarsStencil().size() + interiorDirichletScvfIndexSet_().size();
224
225

        // the local matrices
226
227
228
229
        DynamicMatrix A(numFluxFaces, numFluxFaces, 0.0);
        DynamicMatrix B(numFluxFaces, numPotentials, 0.0);
        DynamicMatrix C(numFaces, numFluxFaces, 0.0);
        DynamicMatrix D(numFaces, numPotentials, 0.0);
230
231
232
233

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

        // solve local system and store matrices
234
        DynamicMatrix copy(B);
235
236
237
238
239
240
241
        A.invert();
        AinvB_ = B.leftmultiply(A);
        CAinv_ = C.rightmultiply(A);
        T_ = multiplyMatrices(CAinv_, copy);
        T_ += D;
    }

242
243
244
245
246
    //! gets local data on a global scvf within the interaction volume
    //! specialization for dim = dimWorld
    template<int d = dim, int dw = dimWorld>
    typename std::enable_if< (d == dw), LocalFaceData >::type
    getLocalFaceData(const SubControlVolumeFace& scvf) const
247
248
249
    {
        auto scvfGlobalIdx = scvf.index();

250
251
        // if interior boundaries are disabled, do simplified search
        if (!enableInteriorBoundaries)
252
        {
253
254
255
256
257
258
259
260
            LocalIndexType localIdx = 0;
            for (auto&& localScvf : localScvfs_)
            {
                if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
                    return LocalFaceData(localIdx, localScvf.insideLocalScvIndex(), false);

                if (!localScvf.boundary() && localScvf.outsideGlobalScvfIndex() == scvfGlobalIdx)
                    return LocalFaceData(localIdx, localScvf.outsideLocalScvIndex(), true);
261

262
263
264
265
266
267
268
269
270
271
272
273
274
                localIdx++;
            }
        }
        else
        {
            // first check the inside data of the faces
            LocalIndexType localIdx = 0;
            for (auto&& localScvf : localScvfs_)
            {
                if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
                    return LocalFaceData(localIdx, localScvf.insideLocalScvIndex(), false);
                localIdx++;
            }
275

276
277
278
279
280
281
282
283
            // then check the outside data of the faces
            localIdx = 0;
            for (auto&& localScvf : localScvfs_)
            {
                if (!localScvf.boundary() && localScvf.outsideGlobalScvfIndex() == scvfGlobalIdx)
                    return LocalFaceData(localIdx, localScvf.outsideLocalScvIndex(), true);
                localIdx++;
            }
284
285
286
287
288
        }

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

289
290
291
292
293
    //! gets local data on a global scvf within the interaction volume
    //! specialization for dim < dimWorld
    template<int d = dim, int dw = dimWorld>
    typename std::enable_if< (d < dw), LocalFaceData >::type
    getLocalFaceData(const SubControlVolumeFace& scvf) const
294
    {
295
        const auto scvfGlobalIdx = scvf.index();
296

297
298
        // if interior boundaries are disabled, do simplified search
        if (!enableInteriorBoundaries)
299
        {
300
301
            LocalIndexType localFaceIdx = 0;
            for (auto&& localScvf : localScvfs_)
302
            {
303
304
                if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
                    return LocalFaceData(localFaceIdx, localScvf.insideLocalScvIndex(), false);
305

306
307
308
309
310
311
312
313
314
315
316
317
318
                if (!localScvf.boundary() && MpfaHelper::contains(localScvf.outsideGlobalScvfIndices(), scvfGlobalIdx))
                {
                    if (localScvf.outsideLocalScvIndices().size() == 1)
                        return LocalFaceData(localFaceIdx, localScvf.outsideLocalScvIndex(), true);

                    // Now we have to find wich local scv is the one the global scvf is embedded in
                    const auto scvGlobalIdx = scvf.insideScvIdx();
                    for (auto localScvIdx : localScvf.outsideLocalScvIndices())
                        if (localScv_(localScvIdx).globalIndex() == scvGlobalIdx)
                            return LocalFaceData(localFaceIdx, localScvIdx, true);

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

320
                localFaceIdx++;
321
            }
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
        }
        else
        {
            // first check the inside data of the faces
            LocalIndexType localFaceIdx = 0;
            for (auto&& localScvf : localScvfs_)
            {
                if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
                    return LocalFaceData(localFaceIdx, localScvf.insideLocalScvIndex(), false);
                localFaceIdx++;
            }

            // then check the outside data
            localFaceIdx = 0;
            for (auto&& localScvf : localScvfs_)
            {
                if (!localScvf.boundary() && MpfaHelper::contains(localScvf.outsideGlobalScvfIndices(), scvfGlobalIdx))
                {
                    if (localScvf.outsideLocalScvIndices().size() == 1)
                        return LocalFaceData(localFaceIdx, localScvf.outsideLocalScvIndex(), true);

                    // Now we have to find wich local scv is the one the global scvf is embedded in
                    const auto scvGlobalIdx = scvf.insideScvIdx();
                    for (auto localScvIdx : localScvf.outsideLocalScvIndices())
                        if (localScv_(localScvIdx).globalIndex() == scvGlobalIdx)
                            return LocalFaceData(localFaceIdx, localScvIdx, true);
348

349
350
351
352
353
                    DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv in the interaction volume for the given scvf with index: " << scvfGlobalIdx);
                }

                localFaceIdx++;
            }
354
355
356
357
358
        }

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

359
    //! Gets the transmissibilities for a sub-control volume face within the interaction volume.
360
361
362
363
364
365
    //! specialization for dim == dimWorld
    template<int d = dim, int dw = dimWorld>
    typename std::enable_if< (d == dw), DynamicVector >::type
    getTransmissibilities(const LocalFaceData& localFaceData) const
    {
        if (localFaceData.isOutside)
366
367
        {
            auto tij = T_[localFaceData.localScvfIndex];
368
            tij *= -1.0;
369
370
371
372
            return tij;
        }
        else
            return T_[localFaceData.localScvfIndex];
373
374
    }

375
376
    //! Gets the transmissibilities for a sub-control volume face within the interaction volume.
    //! specialization for dim < dimWorld.
377
378
379
    template<int d = dim, int dw = dimWorld>
    typename std::enable_if< (d < dw), DynamicVector >::type
    getTransmissibilities(const LocalFaceData& localFaceData) const
380
    {
381
        // If we come from the inside, simply return tij
382
383
384
        if (!localFaceData.isOutside)
            return T_[localFaceData.localScvfIndex];

385
386
        // compute the outside transmissibilities
        DynamicVector tij(volVarsStencil().size() + interiorDirichletScvfIndexSet_().size(), 0.0);
387
388

        // get the local scv and iterate over local coordinates
389
390
391
392
        const auto numLocalScvs = localScvs_.size();
        const auto numDirichletScvfs = dirichletScvfIndexSet_().size();
        const auto& localScv = localScv_(localFaceData.localScvIndex);
        const auto& localScvf = localScvf_(localFaceData.localScvfIndex);
393

394
395
        const auto idxInOutside = this->findIndexInVector(localScvf.outsideLocalScvIndices(), localFaceData.localScvIndex);
        const auto& wijk = wijk_[localFaceData.localScvfIndex][idxInOutside+1];
396
397
        for (int localDir = 0; localDir < dim; localDir++)
        {
398
399
            const auto localScvfIdx = localScv.localScvfIndex(localDir);
            const auto& localScvf = localScvf_(localScvfIdx);
400

401
402
            const auto faceType = localScvf.faceType();
            if (faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet)
403
            {
404
                const auto fluxFaceIndex = this->findIndexInVector(fluxScvfIndexSet_(), localScvfIdx);
405
406
407
408
409
                auto tmp = AinvB_[fluxFaceIndex];
                tmp *= wijk[localDir];

                tij += tmp;
            }
410
            else if (faceType == MpfaFaceTypes::dirichlet)
411
            {
412
                const auto idxInDiriFaces = this->findIndexInVector(dirichletScvfIndexSet_(), localScvfIdx);
413
414
                tij[numLocalScvs + idxInDiriFaces] += wijk[localDir];
            }
415
416
417
418
419
            else if (faceType == MpfaFaceTypes::interiorDirichlet)
            {
                const auto idxInInteriorDiriFaces = this->findIndexInVector(interiorDirichletScvfIndexSet_(), localScvfIdx);
                tij[numLocalScvs + numDirichletScvfs + idxInInteriorDiriFaces] += wijk[localDir];
            }
420
421
422
423
424
425
426
427

            // add entry from the scv unknown
            tij[localFaceData.localScvIndex] -= wijk[localDir];
        }

        return tij;
    }

428
429
430
431
432
433
434
435
436
437
438
    //! Returns the vector of coefficients with which the vector of neumann boundary conditions
    //! has to be multiplied in order to transform them on the scvf this cache belongs to
    DynamicVector getNeumannFluxTransformationCoefficients(const LocalFaceData& localFaceData) const
    {
        auto cij = CAinv_[localFaceData.localScvfIndex];
        if (localFaceData.isOutside)
            cij *= -1.0;
        return cij;
    }

    Scalar getNeumannFlux(const LocalFaceData& localFaceData, unsigned int eqIdx) const
439
    {
440
        if (!onDomainOrInteriorBoundary() || useTpfaBoundary || fluxScvfIndexSet_().size() == 0 )
441
            return 0.0;
442

443
444
445
446
447
448
449
450
451
        // Do the scalar product CAinv_*neumannFluxes[eqIdx]
        assert(CAinv_[localFaceData.localScvfIndex].size() == neumannFluxes_.size() &&
               "Number of columns of matrix does not correspond to number entries in vector!");

        Scalar flux(0.0);
        for (unsigned int i = 0; i < neumannFluxes_.size(); ++i)
            flux += CAinv_[localFaceData.localScvfIndex][i] * neumannFluxes_[i][eqIdx];

        // flip sign if we are coming from the outside
452
        if (localFaceData.isOutside)
453
            return -1.0*flux;
454

455
456
457
        return flux;
    }

458
459
    bool onDomainOrInteriorBoundary() const
    { return onDomainOrInteriorBoundary_; }
460
461
462
463
464
465
466

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

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

467
468
469
    const GlobalIndexSet& globalScvfs() const
    { return globalScvfIndices_; }

470
471
472
    const std::vector<InteriorBoundaryData>& interiorBoundaryData() const
    { return interiorBoundaryData_; }

473
474
private:

475
    const LocalScvfType& localScvf_(const LocalIndexType localScvfIdx) const
476
477
    { return localScvfs_[localScvfIdx]; }

478
    const LocalScvType& localScv_(const LocalIndexType localScvIdx) const
479
480
    { return localScvs_[localScvIdx]; }

481
    const LocalIndexSet& fluxScvfIndexSet_() const
482
483
    { return fluxFaceIndexSet_; }

484
    const LocalIndexSet& dirichletScvfIndexSet_() const
485
486
    { return dirichletFaceIndexSet_; }

487
488
489
490
491
492
    const LocalIndexSet& interiorDirichletScvfIndexSet_() const
    { return interiorDirichletFaceIndexSet_; }

    const LocalIndexSet& interiorBoundaryScvfIndexSet_() const
    { return interiorBoundaryFaceIndexSet_; }

493
    const Element& localElement_(const LocalIndexType localScvIdx) const
494
495
    { return localElements_[localScvIdx]; }

496
    void createLocalEntities_(const Seed& seed)
497
    {
498
499
        const auto numScvs = seed.scvSeeds().size();
        const auto numScvfs = seed.scvfSeeds().size();
500
501
502
503
504

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

Dennis Gläser's avatar
Dennis Gläser committed
505
        for (auto&& scvSeed : seed.scvSeeds())
506
507
        {
            auto element = problem_().model().globalFvGeometry().element(scvSeed.globalIndex());
508
            localScvs_.emplace_back(LocalScvType(problem_(), element, fvGeometry_(), scvSeed));
509
510
511
            localElements_.emplace_back(std::move(element));
        }

Dennis Gläser's avatar
Dennis Gläser committed
512
        for (auto&& scvfSeed : seed.scvfSeeds())
513
514
        {
            // we have to use the "inside" scv face here
Dennis Gläser's avatar
Dennis Gläser committed
515
            auto&& scvf = fvGeometry_().scvf(scvfSeed.insideGlobalScvfIndex());
516
517
518
519
            localScvfs_.emplace_back(LocalScvfType(scvfSeed, scvf));
        }
    }

520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
    void assembleNeumannFluxVector_()
    {
        // initialize the neumann fluxes vector to zero
        neumannFluxes_.resize(fluxFaceIndexSet_.size(), PrimaryVariables(0.0));

        if (!onDomainOrInteriorBoundary() || useTpfaBoundary)
            return;

        LocalIndexType fluxFaceIdx = 0;
        for (auto localFluxFaceIdx : fluxFaceIndexSet_)
        {
            const auto& localScvf = localScvf_(localFluxFaceIdx);
            const auto faceType = localScvf.faceType();

            if (faceType == MpfaFaceTypes::neumann || (faceType == MpfaFaceTypes::interiorNeumann && !facetCoupling))
            {
                const auto& element = localElement_(localScvf.insideLocalScvIndex());
                const auto& globalScvf = fvGeometry_().scvf(localScvf.insideGlobalScvfIndex());
                auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf);
                neumannFlux *= globalScvf.area();
                neumannFlux *= elemVolVars_()[globalScvf.insideScvIdx()].extrusionFactor();

                // The flux is assumed to be prescribed in the form of -D*gradU
                neumannFluxes_[fluxFaceIdx] = neumannFlux;
            }

            fluxFaceIdx++;
        }
    }

550
551
    template<typename GetTensorFunction>
    void assembleLocalMatrices_(const GetTensorFunction& getTensor,
552
553
554
555
                                DynamicMatrix& A,
                                DynamicMatrix& B,
                                DynamicMatrix& C,
                                DynamicMatrix& D)
556
    {
557
558
559
        static const auto xi = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Xi);
        const auto numLocalScvs = localScvs_.size();
        const auto numDirichletScvfs = dirichletScvfIndexSet_().size();
560

561
562
563
        // reserve space for the omegas
        wijk_.resize(localScvfs_.size());

564
        // loop over the local faces
565
        unsigned int rowIdx = 0;
Dennis Gläser's avatar
Dennis Gläser committed
566
        for (auto&& localScvf : localScvfs_)
567
        {
568
569
570
            const auto faceType = localScvf.faceType();
            const bool hasUnknown = faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet;
            const LocalIndexType idxInFluxFaces = hasUnknown ? this->findIndexInVector(fluxScvfIndexSet_(), rowIdx) : -1;
571
572

            // get diffusion tensor in "positive" sub volume
573
574
575
576
577
578
            const auto posLocalScvIdx = localScvf.insideLocalScvIndex();
            const auto& posLocalScv = localScv_(posLocalScvIdx);
            const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
            const auto& posVolVars = elemVolVars_()[posGlobalScv];
            const auto& element = localElement_(posLocalScvIdx);
            const auto tensor = getTensor(element, posVolVars, posGlobalScv);
579
580

            // the omega factors of the "positive" sub volume
581
            auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);
582
            posWijk *= posVolVars.extrusionFactor();
583
584
585
586

            // Check the local directions of the positive sub volume
            for (int localDir = 0; localDir < dim; localDir++)
            {
587
588
589
                const auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir);
                const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
                const auto curFaceType = curLocalScvf.faceType();
590

591
592
                // First, add the entries associated with unknown face pressures
                if (curFaceType != MpfaFaceTypes::dirichlet && curFaceType != MpfaFaceTypes::interiorDirichlet)
593
594
                {
                    // we need the index of the current local scvf in the flux face indices
595
                    auto curIdxInFluxFaces = this->findIndexInVector(fluxScvfIndexSet_(), curLocalScvfIdx);
596

597
                    // this creates an entry in matrix C
598
                    C[rowIdx][curIdxInFluxFaces] += posWijk[localDir];
599
600

                    // proceed depending on if the current face has an unknown
601
                    if (hasUnknown)
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
                    {
                        if (faceType == MpfaFaceTypes::interiorNeumann)
                        {
                            // on interior neumann faces, apply xi factor
                            // However, if this is a boundary face at the same time, don't!
                            A[idxInFluxFaces][curIdxInFluxFaces] += localScvf.globalScvf().boundary() ?
                                                                    posWijk[localDir] : xi*posWijk[localDir];

                            // add values from other domain in case of facet coupling
                            if (facetCoupling && curIdxInFluxFaces == idxInFluxFaces)
                            {
                                // get interior boundary data
                                const auto& data = interiorBoundaryData_[this->findIndexInVector(interiorBoundaryScvfIndexSet_(), curLocalScvfIdx)];

                                // get the volvars on the actual interior boundary face
                                const auto facetVolVars = data.facetVolVars(fvGeometry_());

                                // calculate "leakage factor"
                                const auto n = curLocalScvf.unitOuterNormal();
                                const auto v = [&] ()
                                        {
                                            auto res = n;
                                            res *= -0.5*facetVolVars.extrusionFactor();
                                            res -= curLocalScvf.ip();
                                            res += curLocalScvf.globalScvf().facetCorner();
                                            res /= res.two_norm2();
                                            return res;
                                        } ();

                                // substract from matrix entry
                                A[idxInFluxFaces][curIdxInFluxFaces] -= MpfaHelper::nT_M_v(n, posVolVars.permeability(), v);
                            }
                        }
                        // this means we are on an interior face
                        else
                            A[idxInFluxFaces][curIdxInFluxFaces] += posWijk[localDir];
                    }
639
                }
640
                else if (curFaceType == MpfaFaceTypes::dirichlet)
641
642
                {
                    // the current face is a Dirichlet face and creates entries in D & eventually B
643
                    auto curIdxInDiriFaces = this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx);
644

645
                    D[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
646
647
648
                    if (hasUnknown)
                        B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] -= posWijk[localDir];
                }
649
650
651
652
653
654
655
656
657
                else if (curFaceType == MpfaFaceTypes::interiorDirichlet)
                {
                    // the current face is an interior Dirichlet face and creates entries in D & eventually B
                    auto curIdxInInteriorDiriFaces = this->findIndexInVector(interiorDirichletScvfIndexSet_(), curLocalScvfIdx);

                    D[rowIdx][numLocalScvs + numDirichletScvfs + curIdxInInteriorDiriFaces] += posWijk[localDir];
                    if (hasUnknown)
                        B[idxInFluxFaces][numLocalScvs + numDirichletScvfs + curIdxInInteriorDiriFaces] -= posWijk[localDir];
                }
658
659
660
661
662

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

                if (hasUnknown)
663
664
665
666
667
668
                {
                    if (faceType == MpfaFaceTypes::interiorNeumann && !localScvf.globalScvf().boundary())
                        B[idxInFluxFaces][posLocalScvIdx] += xi*posWijk[localDir];
                    else
                        B[idxInFluxFaces][posLocalScvIdx] += posWijk[localDir];
                }
669
670
            }

671
672
673
            // store the omegas
            wijk_[rowIdx].emplace_back(std::move(posWijk));

674
675
676
            // If we are on an interior or interior neumann face, add values from negative sub volume
            if (!localScvf.globalScvf().boundary() &&
                (faceType == MpfaFaceTypes::interior || faceType == MpfaFaceTypes::interiorNeumann))
677
            {
678
                // loop over all the outside neighbors of this face and add entries
679
                unsigned int indexInOutsideData = 0;
680
                for (auto negLocalScvIdx : localScvf.outsideLocalScvIndices())
681
                {
682
683
684
685
686
                    const auto& negLocalScv = localScv_(negLocalScvIdx);
                    const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
                    const auto& negVolVars = elemVolVars_()[negGlobalScv];
                    const auto& negElement = localElement_(negLocalScvIdx);
                    const auto negTensor = getTensor(negElement, negVolVars, negGlobalScv);
687

688
                    // the omega factors of the "negative" sub volume
689
690
691
692
693
694
                    DimVector negWijk;

                    // if dim < dimWorld, use outside normal vector
                    if (dim < dimWorld)
                    {
                        // outside scvf
695
                        const auto& outsideScvf = fvGeometry_().scvf(localScvf.outsideGlobalScvfIndex(indexInOutsideData));
696
697
698
699
700
701
702
703
                        auto negNormal = outsideScvf.unitOuterNormal();
                        negNormal *= -1.0;
                        negWijk = calculateOmegas_(negLocalScv, negNormal, localScvf.area(), negTensor);
                    }
                    else
                        negWijk = calculateOmegas_(negLocalScv, localScvf.unitOuterNormal(), localScvf.area(), negTensor);

                    // scale by extrusion factpr
704
                    negWijk *= negVolVars.extrusionFactor();
705
706
707

                    // Check local directions of negative sub volume
                    for (int localDir = 0; localDir < dim; localDir++)
708
                    {
709
710
711
                        const auto curLocalScvfIdx = negLocalScv.localScvfIndex(localDir);
                        const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
                        const auto curFaceType = curLocalScvf.faceType();
712

713
                        if (curFaceType != MpfaFaceTypes::dirichlet && curFaceType != MpfaFaceTypes::interiorDirichlet)
714
                        {
715
716
717
718
                            if (faceType == MpfaFaceTypes::interiorNeumann)
                                A[idxInFluxFaces][this->findIndexInVector(fluxScvfIndexSet_(), curLocalScvfIdx)] -= (1-xi)*negWijk[localDir];
                            else
                                A[idxInFluxFaces][this->findIndexInVector(fluxScvfIndexSet_(), curLocalScvfIdx)] -= negWijk[localDir];
719
                        }
720
                        else if (curFaceType == MpfaFaceTypes::interiorDirichlet)
721
                        {
722
723
                            const auto idxInInteriorDiriFaces = this->findIndexInVector(interiorDirichletScvfIndexSet_(), curLocalScvfIdx);
                            B[idxInFluxFaces][numLocalScvs + numDirichletScvfs + idxInInteriorDiriFaces] += negWijk[localDir];
724
                        }
725
726
                        else // Dirichlet face
                            B[idxInFluxFaces][numLocalScvs + this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx)] += negWijk[localDir];
727
728

                        // add entries to matrix B
729
730
731
732
                        if (faceType == MpfaFaceTypes::interiorNeumann)
                            B[idxInFluxFaces][negLocalScvIdx] -= (1-xi)*negWijk[localDir];
                        else
                            B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir];
733
734
                    }

735
736
                    // store the omegas (negative, because of normal vector sign switch)
                    negWijk *= -1.0;
737
                    wijk_[rowIdx].emplace_back(std::move(negWijk));
738
739
740

                    // increment counter in outside data
                    indexInOutsideData++;
741
742
743
744
745
746
747
748
749
750
                }
            }
            // go to the next face
            rowIdx++;
        }
    }

    template<typename GetTensorFunction>
    void assemblePureDirichletSystem_(const GetTensorFunction& getTensor)
    {
751
752
753
754
        const auto numLocalScvs = localScvs_.size();
        const auto numFaces = localScvfs_.size();
        const auto numInteriorDirichletFaces = interiorDirichletScvfIndexSet_().size();
        const auto numPotentials = volVarsStencil().size() + numInteriorDirichletFaces;
755
756
757
758
759
760

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

761
762
763
        // resize the omegas
        wijk_.resize(numFaces);

764
        // Loop over all the faces, in this case these are all dirichlet boundaries
765
        LocalIndexType rowIdx = 0;
Dennis Gläser's avatar
Dennis Gläser committed
766
        for (auto&& localScvf : localScvfs_)
767
768
        {
            // get diffusion tensor in "positive" sub volume
769
770
771
772
773
774
            const auto posLocalScvIdx = localScvf.insideLocalScvIndex();
            const auto& posLocalScv = localScv_(posLocalScvIdx);
            const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
            const auto& posVolVars = elemVolVars_()[posGlobalScv];
            const auto element = localElement_(posLocalScvIdx);
            const auto tensor = getTensor(element, posVolVars, posGlobalScv);
775
776

            // the omega factors of the "positive" sub volume
777
            auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);
778
            posWijk *= posVolVars.extrusionFactor();
779
780
781

            for (int localDir = 0; localDir < dim; localDir++)
            {
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
                // When interior boundaries are disabled, all faces will be of dirichlet type
                if (!enableInteriorBoundaries)
                {
                    const auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir);
                    const auto curIdxInDiriFaces = this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx);
                    T_[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
                    T_[rowIdx][posLocalScvIdx] -= posWijk[localDir];
                }
                else
                {
                    const auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir);
                    const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
                    const auto curFaceType = curLocalScvf.faceType();

                    const auto curIdxInDiriFaces = curFaceType == MpfaFaceTypes::dirichlet ?
                                                   this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx) :
                                                   numInteriorDirichletFaces + this->findIndexInVector(interiorDirichletScvfIndexSet_(), curLocalScvfIdx);
799

800
801
802
                    T_[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
                    T_[rowIdx][posLocalScvIdx] -= posWijk[localDir];
                }
803
804
            }

805
806
807
            // store the omegas
            wijk_[rowIdx].emplace_back(std::move(posWijk));

808
809
810
811
812
            // go to the next face
            rowIdx++;
        }
    }

813
    // TODO: how to do the assertion of positive coefficients for tensors?
814
    DimVector calculateOmegas_(const LocalScvType& localScv,
815
816
                               const GlobalPosition normal,
                               const Scalar area,
817
818
                               const Tensor& T) const
    {
819
        DimVector wijk;
820
        GlobalPosition tmp;
821
822
823
        for (int dir = 0; dir < dim; ++dir)
        {
            T.mv(localScv.innerNormal(dir), tmp);
824
            wijk[dir] = tmp*normal;
825
        }
826
        wijk *= area;
827
        wijk /= localScv.detX();
828
        wijk *= -1.0;
829
830
831
832

        return wijk;
    }

833
    // calculates n_i^T*K_j*nu_k
834
    DimVector calculateOmegas_(const LocalScvType& localScv,
835
836
                               const GlobalPosition normal,
                               const Scalar area,
837
838
                               const Scalar t) const
    {
839
840
841
        // make sure we have positive diffusion coefficients
        assert(t > 0.0 && "non-positive diffusion coefficients cannot be handled by mpfa methods");

842
        DimVector wijk;
843
        GlobalPosition tmp(normal);
844
845
846
847
        tmp *= t;

        for (int dir = 0; dir < dim; ++dir)
            wijk[dir] = tmp*localScv.innerNormal(dir);
848
        wijk *= area;
849
        wijk /= localScv.detX();
850
        wijk *= -1.0;
851
852
853
854
855

        return wijk;
    }

    const Problem& problem_() const
856
    { return *problemPtr_; }
857
858

    const FVElementGeometry& fvGeometry_() const
859
    { return *fvGeometryPtr_; }
860
861

    const ElementVolumeVariables& elemVolVars_() const
862
    { return *elemVolVarsPtr_; }
863

864
865
866
    const Problem* problemPtr_;
    const FVElementGeometry* fvGeometryPtr_;
    const ElementVolumeVariables* elemVolVarsPtr_;
867

868
    bool onDomainOrInteriorBoundary_;
869
870
871
872
873

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

874
875
876
    GlobalIndexSet globalScvfIndices_;
    GlobalIndexSet volVarsStencil_;
    PositionVector volVarsPositions_;
877

878
879
    LocalIndexSet fluxFaceIndexSet_;
    LocalIndexSet dirichletFaceIndexSet_;
880
881
882
    LocalIndexSet interiorDirichletFaceIndexSet_;
    LocalIndexSet interiorBoundaryFaceIndexSet_;
    std::vector<InteriorBoundaryData> interiorBoundaryData_;
883

884
885
    std::vector< std::vector< DimVector > > wijk_;

886
887
888
    DynamicMatrix T_;
    DynamicMatrix AinvB_;
    DynamicMatrix CAinv_;
889

890
    std::vector<PrimaryVariables> neumannFluxes_;
891
892
893
894
895
};

} // end namespace

#endif