interactionvolume.hh 27.6 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
// -*- 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>
27

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

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

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

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

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

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

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

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

65
66
/*!
 * \ingroup Mpfa
67
68
69
70
 * \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.
71
72
 */
template<class TypeTag>
73
74
class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>
          : public CCMpfaOInteractionVolume<TypeTag,
Dennis Gläser's avatar
Dennis Gläser committed
75
                                            CCMpfaOInteractionVolumeTraits<TypeTag>>
76
{
77
    using TraitsType = CCMpfaOInteractionVolumeTraits<TypeTag>;
78
public:
79
80
    // state the traits class type
    using Traits = TraitsType;
81
82
};

Dennis Gläser's avatar
Dennis Gläser committed
83
template<class TypeTag, class Traits>
84
class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase<TypeTag, Traits>
85
{
86
    using ParentType = CCMpfaInteractionVolumeBase<TypeTag, Traits>;
87
88
89
90
91
92
    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 FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
Dennis Gläser's avatar
Dennis Gläser committed
93
    using DualGridNodalIndexSet = typename CCMpfaDualGridIndexSet<TypeTag>::NodalIndexSet;
94

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

Dennis Gläser's avatar
Dennis Gläser committed
101
102
    using Vector = typename Traits::DynamicVector;
    using Matrix = typename Traits::DynamicMatrix;
103
    using Tensor = typename Traits::Tensor;
104

105
106
    using LocalScvType = typename Traits::LocalScvType;
    using LocalScvfType = typename Traits::LocalScvfType;
107

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

115
    using typename ParentType::LocalFaceData;
Dennis Gläser's avatar
Dennis Gläser committed
116
    using typename ParentType::DirichletDataContainer;
117
    using typename ParentType::LocalFaceDataContainer;
Dennis Gläser's avatar
Dennis Gläser committed
118

119
120
121
122
    //! Sets up the local scope for a given iv index set!
    void setUpLocalScope(const IndexSet& indexSet,
                         const Problem& problem,
                         const FVElementGeometry& fvGeometry)
123
    {
124
        //! store a pointer to the index set
Dennis Gläser's avatar
Dennis Gläser committed
125
        indexSetPtr_ = &indexSet;
126

127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
        //! clear previous data
        clear_();

        //! number of interaction-volume-local faces
        numFaces_ = indexSet.numFaces();

        //! number of interaction-volume-local (and node-local) scvs
        const auto& scvIndices = indexSet.globalScvIndices();
        const auto numLocalScvs = indexSet.numScvs();

        //! number of global scvfs appearing in this interaction volume
        const auto numGlobalScvfs = indexSet.nodalIndexSet().numScvfs();

        //! reserve memory for local entities
        elements_.reserve(numLocalScvs);
        scvs_.reserve(numLocalScvs);
        scvfs_.reserve(numFaces_);
        dirichletData_.reserve(numFaces_);
        localFaceData_.reserve(numGlobalScvfs);

        // set up quantities related to sub-control volumes
        for (LocalIndexType scvIdxLocal = 0; scvIdxLocal < numLocalScvs; scvIdxLocal++)
        {
            const auto scvIdxGlobal = scvIndices[scvIdxLocal];
            scvs_.emplace_back(fvGeometry, fvGeometry.scv(scvIdxGlobal), scvIdxLocal, indexSet);
            elements_.emplace_back(fvGeometry.fvGridGeometry().element(scvIdxGlobal));
        }

        // keep track of the number of unknowns etc
        numUnknowns_ = 0;
        numOutsideFaces_ = 0;
        numPotentials_ = numLocalScvs;

        // set up quantitites related to sub-control volume faces
        for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numFaces_; ++faceIdxLocal)
        {
            const auto scvfIdxGlobal = indexSet.scvfIdxGlobal(faceIdxLocal);
            const auto& neighborScvIndicesLocal = indexSet.neighboringLocalScvIndices(faceIdxLocal);
            const auto insideLocalScvIdx = neighborScvIndicesLocal[0];

            // we have to use the "inside" scv face here
            const auto& scvf = fvGeometry.scvf(scvfIdxGlobal);

            // create local face data object for this face
            localFaceData_.emplace_back(faceIdxLocal, insideLocalScvIdx, scvf.index());

            // create iv-local scvf object
            if (scvf.boundary())
            {
                const auto insideElement = elements_[insideLocalScvIdx];
                const auto bcTypes = problem.boundaryTypes(insideElement, scvf);

                if (bcTypes.hasOnlyDirichlet())
                {
                    scvfs_.emplace_back(scvf, neighborScvIndicesLocal, /*isDirichlet*/true, numPotentials_++);
                    dirichletData_.emplace_back(scvf.outsideScvIdx(), scvf.ipGlobal());
                }
                else
                    scvfs_.emplace_back(scvf, neighborScvIndicesLocal, /*isDirichlet*/false, numUnknowns_++);
            }
            else
            {
                scvfs_.emplace_back(scvf, neighborScvIndicesLocal, /*isDirichlet*/false, numUnknowns_++);

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

                    // loop over scvfs in outside scv until we find the one coinciding with current scvf
                    for (int coord = 0; coord < dim; ++coord)
                    {
                        if (indexSet.scvfIdxLocal(outsideLocalScvIdx, coord) == faceIdxLocal)
                        {
                            const auto globalScvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(outsideLocalScvIdx, coord);
                            const auto& flipScvf = fvGeometry.scvf(globalScvfIdx);
                            localFaceData_.emplace_back(faceIdxLocal,         //! iv-local scvf idx
                                                        outsideLocalScvIdx,   //! iv-local scv index
                                                        numOutsideFaces_++,    //! iv-local index in outside faces
                                                        i-1,                  //! scvf-local index in outside faces
                                                        flipScvf.index());   //! global scvf index
                        }
                    }
                }
            }
        }

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

221
222
    //! sets the sizes of the corresponding matrices in the data handle
    void prepareDataHandle(DataHandle& dataHandle)
223
    {
224
225
      // resize the transmissibility matrix in the data handle
      dataHandle.resizeT(numFaces_, numPotentials_);
226

227
228
229
      // resize possible additional containers in the data handle
      if (requireABMatrix_()) dataHandle.resizeAB(numUnknowns_, numPotentials_);
      if (dim < dimWorld) dataHandle.resizeOutsideTij(numOutsideFaces_, numPotentials_);
230
231
    }

Dennis Gläser's avatar
Dennis Gläser committed
232
    //! solves for the transmissibilities subject to a given tensor
233
    template<typename GetTensorFunction>
234
235
236
237
238
    void solveLocalSystem(const GetTensorFunction& getTensor,
                          const Problem& problem,
                          const FVElementGeometry& fvGeometry,
                          const ElementVolumeVariables& elemVolVars,
                          DataHandle& dataHandle)
239
240
    {
        // if only dirichlet faces are present, assemble T_ directly
Dennis Gläser's avatar
Dennis Gläser committed
241
        if (numUnknowns_ == 0)
242
            assemblePureDirichletSystem_(getTensor, problem, fvGeometry, elemVolVars, dataHandle.T());
Dennis Gläser's avatar
Dennis Gläser committed
243
244
245
        else
        {
            // assemble
246
            assembleLocalMatrices_(getTensor, problem, fvGeometry, elemVolVars);
Dennis Gläser's avatar
Dennis Gläser committed
247
248
249
250
251

            // solve
            A_.invert();

            // T = C*A^-1*B + D
252
253
            dataHandle.T() = multiplyMatrices(C_.rightmultiply(A_), B_);
            dataHandle.T() += D_;
Dennis Gläser's avatar
Dennis Gläser committed
254
255
256
257
258
259
260

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

        // // set vol vars stencil & positions pointer in handle
261
        dataHandle.setVolVarsStencilPointer(indexSet().globalScvIndices());
Dennis Gläser's avatar
Dennis Gläser committed
262
263
264
265
266
        dataHandle.setDirichletDataPointer(dirichletData_);

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

Dennis Gläser's avatar
Dennis Gläser committed
269
270
    //! obtain the local data object for a given global scvf
    const LocalFaceData& getLocalFaceData(const SubControlVolumeFace& scvf) const
271
272
273
274
275
276
277
278
279
    {
        //! find corresponding entry in the local face data container
        const auto scvfIdxGlobal = scvf.index();
        auto it = std::find_if(localFaceData_.begin(),
                               localFaceData_.end(),
                               [scvfIdxGlobal] (const LocalFaceData& d) { return d.globalScvfIndex() == scvfIdxGlobal; });
        assert(it != localFaceData_.end() && "Could not find the local face data corresponding to the given scvf");
        return localFaceData_[std::distance(localFaceData_.begin(), it)];
    }
280

281
282
283
    //! returns the grid element corresponding to a given iv-local scv idx
    const Element& element(const LocalIndexType ivLocalScvIdx) const
    { return elements_[ivLocalScvIdx]; }
284

285
286
287
    //! returns the local scvf entity corresponding to a given iv-local scvf idx
    const LocalScvfType& localScvf(const LocalIndexType ivLocalScvfIdx) const
    { return scvfs_[ivLocalScvfIdx]; }
288

289
290
291
    //! returns the local scv entity corresponding to a given iv-local scv idx
    const LocalScvType& localScv(const LocalIndexType ivLocalScvfIdx) const
    { return scvs_[ivLocalScvfIdx]; }
292

293
294
295
    //! returns a reference to the container with the data on Dirichlet boundaries
    const DirichletDataContainer& dirichletData() const
    { return dirichletData_; }
296

297
298
299
    //! returns a reference to the container with the local face data
    const std::vector<LocalFaceData>& localFaceData() const
    { return localFaceData_; }
300

301
    //! returns a reference to the index set of this iv
Dennis Gläser's avatar
Dennis Gläser committed
302
303
    const IndexSet& indexSet() const
    { return *indexSetPtr_; }
304

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

Dennis Gläser's avatar
Dennis Gläser committed
310
311
312
313
314
315
    //! adds the iv index sets living around a vertex to a given container
    //! and stores the the corresponding index in a map for each scvf
    template<class IvIndexSetContainer, class ScvfIndexMap>
    static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
                                              ScvfIndexMap& scvfIndexMap,
                                              const DualGridNodalIndexSet& nodalIndexSet)
316
    {
Dennis Gläser's avatar
Dennis Gläser committed
317
318
        // the global index of the iv index set that is about to be created
        const auto curGlobalIndex = ivIndexSetContainer.size();
319

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

Dennis Gläser's avatar
Dennis Gläser committed
323
324
325
        // store the index mapping
        for (const auto scvfIdx : nodalIndexSet.globalScvfIndices())
            scvfIndexMap[scvfIdx] = curGlobalIndex;
326
327
    }

328
private:
Dennis Gläser's avatar
Dennis Gläser committed
329
330
331
    //! returns a boolean whether or not the AB matrix has to be passed to the handles
    bool requireABMatrix_() const
    {
332
        static const bool requireAB = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Vtk.AddVelocity") || dim < dimWorld;
Dennis Gläser's avatar
Dennis Gläser committed
333
334
        return requireAB;
    }
335

Dennis Gläser's avatar
Dennis Gläser committed
336
    //! clears all the containers
337
    void clear_()
Dennis Gläser's avatar
Dennis Gläser committed
338
339
340
341
    {
        elements_.clear();
        scvs_.clear();
        scvfs_.clear();
342
        localFaceData_.clear();
Dennis Gläser's avatar
Dennis Gläser committed
343
344
        dirichletData_.clear();
    }
345

Dennis Gläser's avatar
Dennis Gläser committed
346
    //! Assembles the local matrices that define the local system of equations and flux expressions
347
    template<typename GetTensorFunction>
348
349
350
351
    void assembleLocalMatrices_(const GetTensorFunction& getTensor,
                                const Problem& problem,
                                const FVElementGeometry& fvGeometry,
                                const ElementVolumeVariables& elemVolVars)
352
    {
Dennis Gläser's avatar
Dennis Gläser committed
353
354
355
356
357
        // reset matrices
        A_ = 0.0;
        B_ = 0.0;
        C_ = 0.0;
        D_ = 0.0;
358

359
        // reserve space for the omegas
Dennis Gläser's avatar
Dennis Gläser committed
360
        wijk_.resize(scvfs_.size());
361

362
        // loop over the local faces
Dennis Gläser's avatar
Dennis Gläser committed
363
        for (unsigned int faceIdx = 0; faceIdx < numFaces_; ++faceIdx)
364
        {
Dennis Gläser's avatar
Dennis Gläser committed
365
            const auto& curLocalScvf = localScvf(faceIdx);
366
            const auto& curGlobalScvf = fvGeometry.scvf(curLocalScvf.globalScvfIndex());
Dennis Gläser's avatar
Dennis Gläser committed
367
368
            const auto curIsDirichlet = curLocalScvf.isDirichlet();
            const auto curLocalDofIdx = curLocalScvf.localDofIndex();
369
370

            // get diffusion tensor in "positive" sub volume
371
            const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
Dennis Gläser's avatar
Dennis Gläser committed
372
373
            const auto posLocalScvIdx = neighborScvIndices[0];
            const auto& posLocalScv = localScv(posLocalScvIdx);
374
375
            const auto& posGlobalScv = fvGeometry.scv(posLocalScv.globalScvIndex());
            const auto& posVolVars = elemVolVars[posGlobalScv];
Dennis Gläser's avatar
Dennis Gläser committed
376
            const auto& posElement = element(posLocalScvIdx);
377
            const auto tensor = getTensor(problem, posElement, posVolVars, fvGeometry, posGlobalScv);
378
379

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

Dennis Gläser's avatar
Dennis Gläser committed
383
384
            // go over the coordinate directions in the positive sub volume
            for (unsigned int localDir = 0; localDir < dim; localDir++)
385
            {
386
                const auto otherLocalScvfIdx = posLocalScv.scvfIdxLocal(localDir);
Dennis Gläser's avatar
Dennis Gläser committed
387
388
                const auto& otherLocalScvf = localScvf(otherLocalScvfIdx);
                const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
389

Dennis Gläser's avatar
Dennis Gläser committed
390
391
392
                // if we are not on a Dirichlet face, add entries associated with unknown face pressures
                // i.e. in matrix C and maybe A (if current face is not a Dirichlet face)
                if (!otherLocalScvf.isDirichlet())
393
                {
Dennis Gläser's avatar
Dennis Gläser committed
394
395
396
                    C_[faceIdx][otherLocalDofIdx] -= posWijk[localDir];
                    if (!curIsDirichlet)
                        A_[curLocalDofIdx][otherLocalDofIdx] -= posWijk[localDir];
397
                }
Dennis Gläser's avatar
Dennis Gläser committed
398
399
                // the current face is a Dirichlet face and creates entries in D & maybe B
                else
400
                {
Dennis Gläser's avatar
Dennis Gläser committed
401
402
403
                    D_[faceIdx][otherLocalDofIdx] -= posWijk[localDir];
                    if (!curIsDirichlet)
                        B_[curLocalDofIdx][otherLocalDofIdx] += posWijk[localDir];
404
                }
405
406

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

Dennis Gläser's avatar
Dennis Gläser committed
410
411
                if (!curIsDirichlet)
                    B_[curLocalDofIdx][posScvLocalDofIdx] -= posWijk[localDir];
412
413
            }

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

Dennis Gläser's avatar
Dennis Gläser committed
417
418
            // If we are on an interior face, add values from negative sub volume
            if (!curGlobalScvf.boundary())
419
            {
420
                // loop over all the outside neighbors of this face and add entries
Dennis Gläser's avatar
Dennis Gläser committed
421
                for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside)
422
                {
Dennis Gläser's avatar
Dennis Gläser committed
423
424
                    const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1];
                    const auto& negLocalScv = localScv(negLocalScvIdx);
425
426
                    const auto& negGlobalScv = fvGeometry.scv(negLocalScv.globalScvIndex());
                    const auto& negVolVars = elemVolVars[negGlobalScv];
Dennis Gläser's avatar
Dennis Gläser committed
427
                    const auto& negElement = element(negLocalScvIdx);
428
                    const auto negTensor = getTensor(problem, negElement, negVolVars, fvGeometry, negGlobalScv);
429

430
                    // the omega factors of the "negative" sub volume
431
432
433
434
435
                    DimVector negWijk;

                    // if dim < dimWorld, use outside normal vector
                    if (dim < dimWorld)
                    {
436
                        const auto& flipScvf = fvGeometry.flipScvf(curGlobalScvf.index(), idxInOutside);
Dennis Gläser's avatar
Dennis Gläser committed
437
                        auto negNormal = flipScvf.unitOuterNormal();
438
                        negNormal *= -1.0;
Dennis Gläser's avatar
Dennis Gläser committed
439
                        negWijk = calculateOmegas_(negLocalScv, negNormal, curGlobalScvf.area(), negTensor);
440
441
                    }
                    else
Dennis Gläser's avatar
Dennis Gläser committed
442
                        negWijk = calculateOmegas_(negLocalScv, curGlobalScvf.unitOuterNormal(), curGlobalScvf.area(), negTensor);
443
444

                    // scale by extrusion factpr
445
                    negWijk *= negVolVars.extrusionFactor();
446

Dennis Gläser's avatar
Dennis Gläser committed
447
                    // go over the coordinate directions in the positive sub volume
448
                    for (int localDir = 0; localDir < dim; localDir++)
449
                    {
450
                        const auto otherLocalScvfIdx = negLocalScv.scvfIdxLocal(localDir);
Dennis Gläser's avatar
Dennis Gläser committed
451
452
                        const auto& otherLocalScvf = localScvf(otherLocalScvfIdx);
                        const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
453

Dennis Gläser's avatar
Dennis Gläser committed
454
455
456
457
                        if (!otherLocalScvf.isDirichlet())
                            A_[curLocalDofIdx][otherLocalDofIdx] += negWijk[localDir];
                        else
                            B_[curLocalDofIdx][otherLocalDofIdx] -= negWijk[localDir];
458
459

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

Dennis Gläser's avatar
Dennis Gläser committed
463
464
                    // store the omegas
                    wijk_[faceIdx].emplace_back(std::move(negWijk));
465
466
467
468
469
                }
            }
        }
    }

Dennis Gläser's avatar
Dennis Gläser committed
470
471
    //! for interaction volumes that have only dirichlet scvfs,
    //! the transmissibility matrix can be assembled directly
472
    template<typename GetTensorFunction>
473
474
475
476
477
    void assemblePureDirichletSystem_(const GetTensorFunction& getTensor,
                                      const Problem& problem,
                                      const FVElementGeometry& fvGeometry,
                                      const ElementVolumeVariables& elemVolVars,
                                      Matrix& T)
478
    {
Dennis Gläser's avatar
Dennis Gläser committed
479
480
        // reset the transmissibility matrix beforehand
        T = 0.0;
481

482
        // Loop over all the faces, in this case these are all dirichlet boundaries
Dennis Gläser's avatar
Dennis Gläser committed
483
        for (unsigned int faceIdx = 0; faceIdx < numFaces_; ++faceIdx)
484
        {
Dennis Gläser's avatar
Dennis Gläser committed
485
            const auto& curLocalScvf = localScvf(faceIdx);
486
            const auto& curGlobalScvf = fvGeometry.scvf(curLocalScvf.globalScvfIndex());
Dennis Gläser's avatar
Dennis Gläser committed
487

488
            // get diffusion tensor in "positive" sub volume
489
            const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
Dennis Gläser's avatar
Dennis Gläser committed
490
491
            const auto posLocalScvIdx = neighborScvIndices[0];
            const auto& posLocalScv = localScv(posLocalScvIdx);
492
493
            const auto& posGlobalScv = fvGeometry.scv(posLocalScv.globalScvIndex());
            const auto& posVolVars = elemVolVars[posGlobalScv];
Dennis Gläser's avatar
Dennis Gläser committed
494
            const auto& posElement = element(posLocalScvIdx);
495
            const auto tensor = getTensor(problem, posElement, posVolVars, fvGeometry, posGlobalScv);
496
497

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

Dennis Gläser's avatar
Dennis Gläser committed
501
            const auto posScvLocalDofIdx = posLocalScv.localDofIndex();
502
503
            for (int localDir = 0; localDir < dim; localDir++)
            {
504
                const auto otherLocalScvfIdx = posLocalScv.scvfIdxLocal(localDir);
Dennis Gläser's avatar
Dennis Gläser committed
505
506
507
508
509
510
511
512
513
514
515
                const auto& otherLocalScvf = localScvf(otherLocalScvfIdx);
                const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
                T[faceIdx][otherLocalDofIdx] -= posWijk[localDir];
                T[faceIdx][posScvLocalDofIdx] += posWijk[localDir];
            }
        }
    }

    //! computes the transmissibilities associated with "outside" faces on surface grids
    void computeOutsideTransmissibilities_(DataHandle& dataHandle) const
    {
516
        assert(dim < dimWorld && "only for dim < dimWorld the outside transmissiblity container has the right size");
Dennis Gläser's avatar
Dennis Gläser committed
517

518
        for (const auto& localFaceData : localFaceData_)
Dennis Gläser's avatar
Dennis Gläser committed
519
        {
520
521
            //! continue only for "outside" faces
            if (!localFaceData.isOutside()) continue;
Dennis Gläser's avatar
Dennis Gläser committed
522

523
524
            const auto localScvIdx = localFaceData.ivLocalInsideScvIndex();
            const auto localScvfIdx = localFaceData.ivLocalScvfIndex();
Dennis Gläser's avatar
Dennis Gläser committed
525
            const auto& posLocalScv = localScv(localScvIdx);
526
            const auto& wijk = wijk_[localScvfIdx][localFaceData.scvfLocalOutsideScvfIndex() + 1];
Dennis Gläser's avatar
Dennis Gläser committed
527

528
529
            //! store the calculated transmissibilities in the data handle
            auto& tij = dataHandle.outsideTij()[localFaceData.ivLocalOutsideScvfIndex()];
Dennis Gläser's avatar
Dennis Gläser committed
530

531
            //! reset transmissibility vector
Dennis Gläser's avatar
Dennis Gläser committed
532
533
            tij = 0.0;

534
            //! add contributions from all local directions
Dennis Gläser's avatar
Dennis Gläser committed
535
536
            for (int localDir = 0; localDir < dim; localDir++)
            {
537
538
539
540
                //! the scvf corresponding to this local direction in the scv
                const auto& curLocalScvf = localScvf(posLocalScv.scvfIdxLocal(localDir));

                //! on interior faces the coefficients of the AB matrix come into play
Dennis Gläser's avatar
Dennis Gläser committed
541
                if (!curLocalScvf.isDirichlet())
542
                {
Dennis Gläser's avatar
Dennis Gläser committed
543
544
545
                    auto tmp = dataHandle.AB()[curLocalScvf.localDofIndex()];
                    tmp *= wijk[localDir];
                    tij -= tmp;
546
547
                }
                else
Dennis Gläser's avatar
Dennis Gläser committed
548
                    tij[curLocalScvf.localDofIndex()] -= wijk[localDir];
549

Dennis Gläser's avatar
Dennis Gläser committed
550
551
                // add entry from the scv unknown
                tij[localScvIdx] += wijk[localDir];
552
553
554
555
            }
        }
    }

Dennis Gläser's avatar
Dennis Gläser committed
556
    // calculates n_i^T*K_j*nu_k
557
    DimVector calculateOmegas_(const LocalScvType& localScv,
558
559
                               const GlobalPosition normal,
                               const Scalar area,
560
561
                               const Tensor& T) const
    {
Dennis Gläser's avatar
Dennis Gläser committed
562
563
564
        // make sure we have positive definite diffsion tensors
        assert(this->tensorIsPositiveDefinite(T) && "only positive definite tensors can be handled by mpfa methods");

565
        DimVector wijk;
566
        GlobalPosition tmp;
567
568
569
        for (int dir = 0; dir < dim; ++dir)
        {
            T.mv(localScv.innerNormal(dir), tmp);
570
            wijk[dir] = tmp*normal;
571
        }
572
        wijk *= area;
573
574
575
576
577
        wijk /= localScv.detX();

        return wijk;
    }

578
    // calculates n_i^T*K_j*nu_k
579
    DimVector calculateOmegas_(const LocalScvType& localScv,
580
581
                               const GlobalPosition normal,
                               const Scalar area,
582
583
                               const Scalar t) const
    {
584
585
586
        // make sure we have positive diffusion coefficients
        assert(t > 0.0 && "non-positive diffusion coefficients cannot be handled by mpfa methods");

587
        DimVector wijk;
588
        GlobalPosition tmp(normal);
589
590
591
592
        tmp *= t;

        for (int dir = 0; dir < dim; ++dir)
            wijk[dir] = tmp*localScv.innerNormal(dir);
593
        wijk *= area;
594
595
596
597
598
        wijk /= localScv.detX();

        return wijk;
    }

Dennis Gläser's avatar
Dennis Gläser committed
599
    const IndexSet* indexSetPtr_;
600

601
    // Variables defining the local scope
Dennis Gläser's avatar
Dennis Gläser committed
602
603
604
    std::vector<Element> elements_;
    std::vector<LocalScvType> scvs_;
    std::vector<LocalScvfType> scvfs_;
605
    std::vector<LocalFaceData> localFaceData_;
Dennis Gläser's avatar
Dennis Gläser committed
606
    DirichletDataContainer dirichletData_;
607

Dennis Gläser's avatar
Dennis Gläser committed
608
609
610
611
    // sizes involved in the local matrices
    unsigned int numFaces_;
    unsigned int numUnknowns_;
    unsigned int numPotentials_;
612
    unsigned int numOutsideFaces_;
613

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

Dennis Gläser's avatar
Dennis Gläser committed
619
620
621
622
623
    // Matrices involved in transmissibility calculations
    Matrix A_;
    Matrix B_;
    Matrix C_;
    Matrix D_;
624
625
626
627
628
};

} // end namespace

#endif