interactionvolume.hh 37.4 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
    using typename ParentType::GlobalLocalFaceDataPair;
132
133
    using typename ParentType::LocalIndexType;
    using typename ParentType::LocalIndexSet;
134
    using typename ParentType::LocalFaceData;
135
136
137
138
139
    using typename ParentType::GlobalIndexSet;
    using typename ParentType::PositionVector;
    using typename ParentType::Seed;

    CCMpfaOInteractionVolume(const Seed& seed,
140
141
142
                             const Problem& problem,
                             const FVElementGeometry& fvGeometry,
                             const ElementVolumeVariables& elemVolVars)
143
144
145
    : problemPtr_(&problem),
      fvGeometryPtr_(&fvGeometry),
      elemVolVarsPtr_(&elemVolVars),
146
      onDomainOrInteriorBoundary_(seed.onDomainOrInteriorBoundary())
147
    {
148
149
        // set up the local scope of this interaction volume
        bind(seed);
150

151
152
        // initialize the vector containing the neumann fluxes
        assembleNeumannFluxVector_();
153
154
155
156
157
    }

    template<typename GetTensorFunction>
    void solveLocalSystem(const GetTensorFunction& getTensor)
    {
158
        const auto numFluxFaces = fluxScvfIndexSet_().size();
159
160
161
162
163

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

164
165
        const auto numFaces = localScvfs_.size();
        const auto numPotentials = volVarsStencil().size() + interiorDirichletScvfIndexSet_().size();
166
167

        // the local matrices
168
169
170
171
        DynamicMatrix A(numFluxFaces, numFluxFaces, 0.0);
        DynamicMatrix B(numFluxFaces, numPotentials, 0.0);
        DynamicMatrix C(numFaces, numFluxFaces, 0.0);
        DynamicMatrix D(numFaces, numPotentials, 0.0);
172
173
174
175

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

        // solve local system and store matrices
176
        DynamicMatrix copy(B);
177
178
179
180
181
182
183
        A.invert();
        AinvB_ = B.leftmultiply(A);
        CAinv_ = C.rightmultiply(A);
        T_ = multiplyMatrices(CAinv_, copy);
        T_ += D;
    }

184
    //! Gets the transmissibilities for a sub-control volume face within the interaction volume.
185
186
187
188
189
190
    //! 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)
191
192
        {
            auto tij = T_[localFaceData.localScvfIndex];
193
            tij *= -1.0;
194
195
196
197
            return tij;
        }
        else
            return T_[localFaceData.localScvfIndex];
198
199
    }

200
201
    //! Gets the transmissibilities for a sub-control volume face within the interaction volume.
    //! specialization for dim < dimWorld.
202
203
204
    template<int d = dim, int dw = dimWorld>
    typename std::enable_if< (d < dw), DynamicVector >::type
    getTransmissibilities(const LocalFaceData& localFaceData) const
205
    {
206
        // If we come from the inside, simply return tij
207
208
209
        if (!localFaceData.isOutside)
            return T_[localFaceData.localScvfIndex];

210
211
        // compute the outside transmissibilities
        DynamicVector tij(volVarsStencil().size() + interiorDirichletScvfIndexSet_().size(), 0.0);
212
213

        // get the local scv and iterate over local coordinates
214
215
216
217
        const auto numLocalScvs = localScvs_.size();
        const auto numDirichletScvfs = dirichletScvfIndexSet_().size();
        const auto& localScv = localScv_(localFaceData.localScvIndex);
        const auto& localScvf = localScvf_(localFaceData.localScvfIndex);
218

219
220
        const auto idxInOutside = this->findIndexInVector(localScvf.outsideLocalScvIndices(), localFaceData.localScvIndex);
        const auto& wijk = wijk_[localFaceData.localScvfIndex][idxInOutside+1];
221
222
        for (int localDir = 0; localDir < dim; localDir++)
        {
223
224
            const auto localScvfIdx = localScv.localScvfIndex(localDir);
            const auto& localScvf = localScvf_(localScvfIdx);
225

226
227
            const auto faceType = localScvf.faceType();
            if (faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet)
228
            {
229
                const auto fluxFaceIndex = this->findIndexInVector(fluxScvfIndexSet_(), localScvfIdx);
230
231
232
233
234
                auto tmp = AinvB_[fluxFaceIndex];
                tmp *= wijk[localDir];

                tij += tmp;
            }
235
            else if (faceType == MpfaFaceTypes::dirichlet)
236
            {
237
                const auto idxInDiriFaces = this->findIndexInVector(dirichletScvfIndexSet_(), localScvfIdx);
238
239
                tij[numLocalScvs + idxInDiriFaces] += wijk[localDir];
            }
240
241
242
243
244
            else if (faceType == MpfaFaceTypes::interiorDirichlet)
            {
                const auto idxInInteriorDiriFaces = this->findIndexInVector(interiorDirichletScvfIndexSet_(), localScvfIdx);
                tij[numLocalScvs + numDirichletScvfs + idxInInteriorDiriFaces] += wijk[localDir];
            }
245
246
247
248
249
250
251
252

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

        return tij;
    }

253
254
255
256
257
258
259
260
261
262
263
    //! 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
264
    {
265
        if (!onDomainOrInteriorBoundary() || useTpfaBoundary || fluxScvfIndexSet_().size() == 0 )
266
            return 0.0;
267

268
269
270
271
272
273
274
275
276
        // 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
277
        if (localFaceData.isOutside)
278
            return -1.0*flux;
279

280
281
282
        return flux;
    }

283
284
285
    const LocalFaceData& getLocalFaceData(const SubControlVolumeFace& scvf) const
    { return globalLocalScvfPairedData_[this->findIndexInVector(globalScvfIndices_, scvf.index())].second; }

286
287
    bool onDomainOrInteriorBoundary() const
    { return onDomainOrInteriorBoundary_; }
288
289
290
291
292
293
294

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

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

295
296
    const std::vector<GlobalLocalFaceDataPair>& globalLocalScvfPairedData() const
    { return globalLocalScvfPairedData_; }
297

298
299
300
    const std::vector<InteriorBoundaryData>& interiorBoundaryData() const
    { return interiorBoundaryData_; }

301
302
private:

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

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

309
    const LocalIndexSet& fluxScvfIndexSet_() const
310
311
    { return fluxFaceIndexSet_; }

312
    const LocalIndexSet& dirichletScvfIndexSet_() const
313
314
    { return dirichletFaceIndexSet_; }

315
316
317
318
319
320
    const LocalIndexSet& interiorDirichletScvfIndexSet_() const
    { return interiorDirichletFaceIndexSet_; }

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

321
    const Element& localElement_(const LocalIndexType localScvIdx) const
322
323
    { return localElements_[localScvIdx]; }

324
    void bind(const Seed& seed)
325
    {
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
        const auto numLocalScvs = seed.scvSeeds().size();
        const auto numLocalScvfs = seed.scvfSeeds().size();
        const auto numGlobalScvfs = seed.globalScvfIndices().size();
        const auto maxNumVolVars = numLocalScvs + numLocalScvfs;

        //! reserve memory for local entities
        localElements_.reserve(numLocalScvs);
        localScvs_.reserve(numLocalScvs);
        localScvfs_.reserve(numLocalScvfs);
        globalLocalScvfPairedData_.reserve(numGlobalScvfs);
        globalScvfIndices_.reserve(numGlobalScvfs);

        //! reserve memory for the index sets
        volVarsStencil_ = seed.globalScvIndices(); // boundary vol vars are placed at the end
        volVarsStencil_.reserve(maxNumVolVars);
        volVarsPositions_.reserve(maxNumVolVars);
        dirichletFaceIndexSet_.reserve(numLocalScvfs);
        interiorDirichletFaceIndexSet_.reserve(numLocalScvfs);
        interiorBoundaryFaceIndexSet_.reserve(numLocalScvfs);
        interiorBoundaryData_.reserve(numLocalScvfs);
        fluxFaceIndexSet_.reserve(numLocalScvfs);
347

348
        // set up quantities related to sub-control volumes
Dennis Gläser's avatar
Dennis Gläser committed
349
        for (auto&& scvSeed : seed.scvSeeds())
350
        {
351
352
            const auto element = problem_().model().globalFvGeometry().element(scvSeed.globalIndex());
            localScvs_.emplace_back(problem_(), element, fvGeometry_(), scvSeed);
353
            localElements_.emplace_back(std::move(element));
354
            volVarsPositions_.push_back(localScvs_.back().center());
355
356
        }

357
358
        // set up quantitites related to sub-control volume faces
        LocalIndexType localFaceIdx = 0;
Dennis Gläser's avatar
Dennis Gläser committed
359
        for (auto&& scvfSeed : seed.scvfSeeds())
360
        {
361
362
            const auto faceType = scvfSeed.faceType();

363
            // we have to use the "inside" scv face here
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
            const auto& scvf = fvGeometry_().scvf(scvfSeed.insideGlobalScvfIndex());
            localScvfs_.emplace_back(scvfSeed, scvf);

            // create global/local face data for this face
            // we simultaneously store the corresponding global scvf indices (allows global to local mapping later)
            globalLocalScvfPairedData_.emplace_back(&scvf, LocalFaceData(localFaceIdx, scvfSeed.insideLocalScvIndex(), false));
            globalScvfIndices_.push_back(scvf.index());

            // set data depending on the face type
            // obtain the local scvf entity just inserted
            const auto& localScvf = localScvfs_.back();

            // interior faces are flux faces
            // also, set up outside global/local data for all the neighbors
            if (faceType == MpfaFaceTypes::interior)
            {
                for (unsigned int i = 0; i < localScvf.outsideGlobalScvfIndices().size(); ++i)
                {
                    const auto& outsideScvf = fvGeometry_().scvf(localScvf.outsideGlobalScvfIndex(i));
                    globalLocalScvfPairedData_.emplace_back( &outsideScvf,
                                                             LocalFaceData(localFaceIdx, scvfSeed.outsideLocalScvIndex(i), true) );
                    globalScvfIndices_.push_back(outsideScvf.index());
                }
                fluxFaceIndexSet_.push_back(localFaceIdx++);
            }
            // dirichlet faces are in the "stencil"
            else if (faceType == MpfaFaceTypes::dirichlet)
            {
                volVarsStencil_.push_back(localScvf.outsideGlobalScvIndex());
                volVarsPositions_.push_back(localScvf.ip());
                dirichletFaceIndexSet_.push_back(localFaceIdx++);
            }
            // neumann faces have an unknown associated with them
            else if (faceType == MpfaFaceTypes::neumann)
            {
                fluxFaceIndexSet_.push_back(localFaceIdx++);
            }
            // interior neumann faces additionally produce interior boundary data
            else if (faceType == MpfaFaceTypes::interiorNeumann)
            {
                interiorBoundaryData_.push_back(InteriorBoundaryData(problem_(),
                                                                     localScvf.insideGlobalScvIndex(),
                                                                     localScvf.insideGlobalScvfIndex(),
                                                                     fluxFaceIndexSet_.size(),
                                                                     faceType));
                fluxFaceIndexSet_.push_back(localFaceIdx);
                interiorBoundaryFaceIndexSet_.push_back(localFaceIdx++);
            }
            // as well as interior Dirichlet boundaries
            else if (faceType == MpfaFaceTypes::interiorDirichlet)
            {
                interiorBoundaryData_.push_back(InteriorBoundaryData(problem_(),
                                                                     localScvf.insideGlobalScvIndex(),
                                                                     localScvf.insideGlobalScvfIndex(),
                                                                     interiorDirichletFaceIndexSet_.size(),
                                                                     faceType));
                interiorDirichletFaceIndexSet_.push_back(localFaceIdx);
                interiorBoundaryFaceIndexSet_.push_back(localFaceIdx++);
            }
            else
                DUNE_THROW(Dune::InvalidStateException, "Face type can not be handled by the mpfa o-method.");
425
426
427
        }
    }

428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
    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++;
        }
    }

458
459
    template<typename GetTensorFunction>
    void assembleLocalMatrices_(const GetTensorFunction& getTensor,
460
461
462
463
                                DynamicMatrix& A,
                                DynamicMatrix& B,
                                DynamicMatrix& C,
                                DynamicMatrix& D)
464
    {
465
466
467
        static const auto xi = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Xi);
        const auto numLocalScvs = localScvs_.size();
        const auto numDirichletScvfs = dirichletScvfIndexSet_().size();
468

469
470
471
        // reserve space for the omegas
        wijk_.resize(localScvfs_.size());

472
        // loop over the local faces
473
        unsigned int rowIdx = 0;
Dennis Gläser's avatar
Dennis Gläser committed
474
        for (auto&& localScvf : localScvfs_)
475
        {
476
477
478
            const auto faceType = localScvf.faceType();
            const bool hasUnknown = faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet;
            const LocalIndexType idxInFluxFaces = hasUnknown ? this->findIndexInVector(fluxScvfIndexSet_(), rowIdx) : -1;
479
480

            // get diffusion tensor in "positive" sub volume
481
482
483
484
485
            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);
486
            const auto tensor = getTensor(problem_(), element, posVolVars, fvGeometry_(), posGlobalScv);
487
488

            // the omega factors of the "positive" sub volume
489
            auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);
490
            posWijk *= posVolVars.extrusionFactor();
491
492
493
494

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

499
500
                // First, add the entries associated with unknown face pressures
                if (curFaceType != MpfaFaceTypes::dirichlet && curFaceType != MpfaFaceTypes::interiorDirichlet)
501
502
                {
                    // we need the index of the current local scvf in the flux face indices
503
                    auto curIdxInFluxFaces = this->findIndexInVector(fluxScvfIndexSet_(), curLocalScvfIdx);
504

505
                    // this creates an entry in matrix C
506
                    C[rowIdx][curIdxInFluxFaces] += posWijk[localDir];
507
508

                    // proceed depending on if the current face has an unknown
509
                    if (hasUnknown)
510
511
512
513
514
515
516
517
518
519
520
521
522
523
                    {
                        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)];

524
525
                                // obtain the complete data on the facet element
                                const auto completeFacetData = data.completeCoupledFacetData();
526
527
528
529
530
531

                                // calculate "leakage factor"
                                const auto n = curLocalScvf.unitOuterNormal();
                                const auto v = [&] ()
                                        {
                                            auto res = n;
532
                                            res *= -0.5*completeFacetData.volVars().extrusionFactor();
533
534
535
536
537
538
                                            res -= curLocalScvf.ip();
                                            res += curLocalScvf.globalScvf().facetCorner();
                                            res /= res.two_norm2();
                                            return res;
                                        } ();

539
540
541
542
543
544
545
                                // substract n*T*v from diagonal matrix entry
                                const auto facetTensor = getTensor(completeFacetData.problem(),
                                                                   completeFacetData.element(),
                                                                   completeFacetData.volVars(),
                                                                   completeFacetData.fvGeometry(),
                                                                   completeFacetData.scv());
                                A[idxInFluxFaces][curIdxInFluxFaces] -= MpfaHelper::nT_M_v(n, facetTensor, v);
546
547
548
549
550
551
                            }
                        }
                        // this means we are on an interior face
                        else
                            A[idxInFluxFaces][curIdxInFluxFaces] += posWijk[localDir];
                    }
552
                }
553
                else if (curFaceType == MpfaFaceTypes::dirichlet)
554
555
                {
                    // the current face is a Dirichlet face and creates entries in D & eventually B
556
                    auto curIdxInDiriFaces = this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx);
557

558
                    D[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
559
560
561
                    if (hasUnknown)
                        B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] -= posWijk[localDir];
                }
562
563
564
565
566
567
568
569
570
                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];
                }
571
572
573
574
575

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

                if (hasUnknown)
576
577
578
579
580
581
                {
                    if (faceType == MpfaFaceTypes::interiorNeumann && !localScvf.globalScvf().boundary())
                        B[idxInFluxFaces][posLocalScvIdx] += xi*posWijk[localDir];
                    else
                        B[idxInFluxFaces][posLocalScvIdx] += posWijk[localDir];
                }
582
583
            }

584
585
586
            // store the omegas
            wijk_[rowIdx].emplace_back(std::move(posWijk));

587
588
589
            // 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))
590
            {
591
                // loop over all the outside neighbors of this face and add entries
592
                unsigned int indexInOutsideData = 0;
593
                for (auto negLocalScvIdx : localScvf.outsideLocalScvIndices())
594
                {
595
596
597
598
                    const auto& negLocalScv = localScv_(negLocalScvIdx);
                    const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
                    const auto& negVolVars = elemVolVars_()[negGlobalScv];
                    const auto& negElement = localElement_(negLocalScvIdx);
599
                    const auto negTensor = getTensor(problem_(), negElement, negVolVars, fvGeometry_(), negGlobalScv);
600

601
                    // the omega factors of the "negative" sub volume
602
603
604
605
606
607
                    DimVector negWijk;

                    // if dim < dimWorld, use outside normal vector
                    if (dim < dimWorld)
                    {
                        // outside scvf
608
                        const auto& outsideScvf = fvGeometry_().scvf(localScvf.outsideGlobalScvfIndex(indexInOutsideData));
609
610
611
612
613
614
615
616
                        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
617
                    negWijk *= negVolVars.extrusionFactor();
618
619
620

                    // Check local directions of negative sub volume
                    for (int localDir = 0; localDir < dim; localDir++)
621
                    {
622
623
624
                        const auto curLocalScvfIdx = negLocalScv.localScvfIndex(localDir);
                        const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
                        const auto curFaceType = curLocalScvf.faceType();
625

626
                        if (curFaceType != MpfaFaceTypes::dirichlet && curFaceType != MpfaFaceTypes::interiorDirichlet)
627
                        {
628
629
630
631
                            if (faceType == MpfaFaceTypes::interiorNeumann)
                                A[idxInFluxFaces][this->findIndexInVector(fluxScvfIndexSet_(), curLocalScvfIdx)] -= (1-xi)*negWijk[localDir];
                            else
                                A[idxInFluxFaces][this->findIndexInVector(fluxScvfIndexSet_(), curLocalScvfIdx)] -= negWijk[localDir];
632
                        }
633
                        else if (curFaceType == MpfaFaceTypes::interiorDirichlet)
634
                        {
635
636
                            const auto idxInInteriorDiriFaces = this->findIndexInVector(interiorDirichletScvfIndexSet_(), curLocalScvfIdx);
                            B[idxInFluxFaces][numLocalScvs + numDirichletScvfs + idxInInteriorDiriFaces] += negWijk[localDir];
637
                        }
638
639
                        else // Dirichlet face
                            B[idxInFluxFaces][numLocalScvs + this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx)] += negWijk[localDir];
640
641

                        // add entries to matrix B
642
643
644
645
                        if (faceType == MpfaFaceTypes::interiorNeumann)
                            B[idxInFluxFaces][negLocalScvIdx] -= (1-xi)*negWijk[localDir];
                        else
                            B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir];
646
647
                    }

648
649
                    // store the omegas (negative, because of normal vector sign switch)
                    negWijk *= -1.0;
650
                    wijk_[rowIdx].emplace_back(std::move(negWijk));
651
652
653

                    // increment counter in outside data
                    indexInOutsideData++;
654
655
656
657
658
659
660
661
662
663
                }
            }
            // go to the next face
            rowIdx++;
        }
    }

    template<typename GetTensorFunction>
    void assemblePureDirichletSystem_(const GetTensorFunction& getTensor)
    {
664
665
666
667
        const auto numLocalScvs = localScvs_.size();
        const auto numFaces = localScvfs_.size();
        const auto numInteriorDirichletFaces = interiorDirichletScvfIndexSet_().size();
        const auto numPotentials = volVarsStencil().size() + numInteriorDirichletFaces;
668
669
670
671
672
673

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

674
675
676
        // resize the omegas
        wijk_.resize(numFaces);

677
        // Loop over all the faces, in this case these are all dirichlet boundaries
678
        LocalIndexType rowIdx = 0;
Dennis Gläser's avatar
Dennis Gläser committed
679
        for (auto&& localScvf : localScvfs_)
680
681
        {
            // get diffusion tensor in "positive" sub volume
682
683
684
685
686
            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);
687
            const auto tensor = getTensor(problem_(), element, posVolVars, fvGeometry_(), posGlobalScv);
688
689

            // the omega factors of the "positive" sub volume
690
            auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);
691
            posWijk *= posVolVars.extrusionFactor();
692
693
694

            for (int localDir = 0; localDir < dim; localDir++)
            {
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
                // 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);
712

713
714
715
                    T_[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
                    T_[rowIdx][posLocalScvIdx] -= posWijk[localDir];
                }
716
717
            }

718
719
720
            // store the omegas
            wijk_[rowIdx].emplace_back(std::move(posWijk));

721
722
723
724
725
            // go to the next face
            rowIdx++;
        }
    }

726
    // TODO: how to do the assertion of positive coefficients for tensors?
727
    DimVector calculateOmegas_(const LocalScvType& localScv,
728
729
                               const GlobalPosition normal,
                               const Scalar area,
730
731
                               const Tensor& T) const
    {
732
        DimVector wijk;
733
        GlobalPosition tmp;
734
735
736
        for (int dir = 0; dir < dim; ++dir)
        {
            T.mv(localScv.innerNormal(dir), tmp);
737
            wijk[dir] = tmp*normal;
738
        }
739
        wijk *= area;
740
        wijk /= localScv.detX();
741
        wijk *= -1.0;
742
743
744
745

        return wijk;
    }

746
    // calculates n_i^T*K_j*nu_k
747
    DimVector calculateOmegas_(const LocalScvType& localScv,
748
749
                               const GlobalPosition normal,
                               const Scalar area,
750
751
                               const Scalar t) const
    {
752
753
754
        // make sure we have positive diffusion coefficients
        assert(t > 0.0 && "non-positive diffusion coefficients cannot be handled by mpfa methods");

755
        DimVector wijk;
756
        GlobalPosition tmp(normal);
757
758
759
760
        tmp *= t;

        for (int dir = 0; dir < dim; ++dir)
            wijk[dir] = tmp*localScv.innerNormal(dir);
761
        wijk *= area;
762
        wijk /= localScv.detX();
763
        wijk *= -1.0;
764
765
766
767
768

        return wijk;
    }

    const Problem& problem_() const
769
    { return *problemPtr_; }
770
771

    const FVElementGeometry& fvGeometry_() const
772
    { return *fvGeometryPtr_; }
773
774

    const ElementVolumeVariables& elemVolVars_() const
775
    { return *elemVolVarsPtr_; }
776

777
778
779
    const Problem* problemPtr_;
    const FVElementGeometry* fvGeometryPtr_;
    const ElementVolumeVariables* elemVolVarsPtr_;
780

781
    bool onDomainOrInteriorBoundary_;
782

783
    // Variables defining the local scope
784
785
786
    std::vector<Element> localElements_;
    std::vector<LocalScvType> localScvs_;
    std::vector<LocalScvfType> localScvfs_;
787
    std::vector<GlobalLocalFaceDataPair> globalLocalScvfPairedData_;
788
    GlobalIndexSet globalScvfIndices_;
789

790
791
    GlobalIndexSet volVarsStencil_;
    PositionVector volVarsPositions_;
792

793
794
    LocalIndexSet fluxFaceIndexSet_;
    LocalIndexSet dirichletFaceIndexSet_;
795
796
797
    LocalIndexSet interiorDirichletFaceIndexSet_;
    LocalIndexSet interiorBoundaryFaceIndexSet_;
    std::vector<InteriorBoundaryData> interiorBoundaryData_;
798

799
    // Quantities depending on the tensor the system is solved for
800
    std::vector< std::vector< DimVector > > wijk_;
801
802
803
    DynamicMatrix T_;
    DynamicMatrix AinvB_;
    DynamicMatrix CAinv_;
804

805
    // stores all the neumann fluxes appearing in this interaction volume
806
    std::vector<PrimaryVariables> neumannFluxes_;
807
808
809
810
811
};

} // end namespace

#endif