| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup CCMpfaDiscretization | ||
| 10 | * \brief Class for the assembly of the local systems of equations | ||
| 11 | * involved in the transmissibility computation in an mpfa-o type manner. | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH | ||
| 14 | #define DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH | ||
| 15 | |||
| 16 | #include <algorithm> | ||
| 17 | |||
| 18 | #include <dumux/discretization/cellcentered/mpfa/localassemblerbase.hh> | ||
| 19 | #include <dumux/discretization/cellcentered/mpfa/localassemblerhelper.hh> | ||
| 20 | #include <dumux/discretization/cellcentered/mpfa/computetransmissibility.hh> | ||
| 21 | |||
| 22 | namespace Dumux { | ||
| 23 | |||
| 24 | /*! | ||
| 25 | * \ingroup CCMpfaDiscretization | ||
| 26 | * \brief Specialization of the interaction volume-local | ||
| 27 | * assembler class for the schemes using an mpfa-o type assembly. | ||
| 28 | * | ||
| 29 | * \tparam P The problem type | ||
| 30 | * \tparam EG The element finite volume geometry | ||
| 31 | * \tparam EV The element volume variables type | ||
| 32 | */ | ||
| 33 | template< class P, class EG, class EV > | ||
| 34 | class MpfaOInteractionVolumeAssembler | ||
| 35 | : public InteractionVolumeAssemblerBase< P, EG, EV > | ||
| 36 | { | ||
| 37 | using ParentType = InteractionVolumeAssemblerBase< P, EG, EV >; | ||
| 38 | using Helper = InteractionVolumeAssemblerHelper; | ||
| 39 | |||
| 40 | template< class IV > | ||
| 41 | using Scalar = typename IV::Traits::MatVecTraits::FaceVector::value_type; | ||
| 42 | |||
| 43 | public: | ||
| 44 | //! Pull up constructor of the base class | ||
| 45 |
10/10✓ Branch 0 taken 398494 times.
✓ Branch 1 taken 47445356 times.
✓ Branch 2 taken 398498 times.
✓ Branch 3 taken 48192552 times.
✓ Branch 4 taken 2404514 times.
✓ Branch 5 taken 12812395 times.
✓ Branch 6 taken 2404512 times.
✓ Branch 7 taken 12118861 times.
✓ Branch 8 taken 2 times.
✓ Branch 9 taken 53662 times.
|
134823144 | using ParentType::ParentType; |
| 46 | |||
| 47 | /*! | ||
| 48 | * \brief Assembles the matrices involved in the flux | ||
| 49 | * expressions and the local system of equations | ||
| 50 | * within an interaction volume in an mpfa-o type way. | ||
| 51 | * | ||
| 52 | * \param handle The data handle in which the matrices are stored | ||
| 53 | * \param iv The interaction volume | ||
| 54 | * \param getT Lambda to evaluate the scv-wise tensors | ||
| 55 | * \param wijZeroThresh Threshold below which transmissibilities are taken to be zero. | ||
| 56 | * On the basis of this threshold, trivial (0 = 0) rows in the | ||
| 57 | * A matrix are identified and modified accordingly in order to | ||
| 58 | * avoid ending up with singular matrices. This can occur when the | ||
| 59 | * tensor is zero in some cells. | ||
| 60 | */ | ||
| 61 | template< class DataHandle, class IV, class TensorFunc > | ||
| 62 | 78449942 | void assembleMatrices(DataHandle& handle, IV& iv, const TensorFunc& getT, Scalar<IV> wijZeroThresh = 0.0) | |
| 63 | { | ||
| 64 |
2/4✓ Branch 6 taken 19515906 times.
✓ Branch 7 taken 24208443 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
549149594 | const auto zeroRows = assembleLocalMatrices_(handle.A(), handle.AB(), handle.CA(), handle.T(), handle.omegas(), iv, getT, wijZeroThresh); |
| 65 | |||
| 66 | // maybe solve the local system | ||
| 67 |
2/2✓ Branch 0 taken 43017228 times.
✓ Branch 1 taken 192 times.
|
77036084 | if (iv.numUnknowns() > 0) |
| 68 | { | ||
| 69 | // for now, B is carried in what will be AB after solve | ||
| 70 | 78449750 | auto& B = handle.AB(); | |
| 71 | 78449750 | auto& A = handle.A(); | |
| 72 | |||
| 73 | // If the tensor is zero in some cells, we might have zero rows in the matrix. | ||
| 74 | // In this case we set the diagonal entry of A to 1.0 and the row of B to zero. | ||
| 75 | // This will ultimatively lead to zero transmissibilities for this face. | ||
| 76 |
4/4✓ Branch 0 taken 75879569 times.
✓ Branch 1 taken 43724157 times.
✓ Branch 2 taken 75879569 times.
✓ Branch 3 taken 43724157 times.
|
499599326 | for (const auto& zeroRowIndices : zeroRows) |
| 77 | { | ||
| 78 | 132125038 | const auto zeroRowDofIdx = zeroRowIndices.first; | |
| 79 |
2/2✓ Branch 0 taken 75879569 times.
✓ Branch 1 taken 301423121 times.
|
792935854 | for (auto& row : A) |
| 80 | 1057371556 | row[zeroRowDofIdx] = 0.0; | |
| 81 | 264250076 | A[zeroRowDofIdx] = 0.0; | |
| 82 | 264250076 | A[zeroRowDofIdx][zeroRowDofIdx] = 1.0; | |
| 83 | 396375114 | B[zeroRowDofIdx] = 0.0; | |
| 84 | } | ||
| 85 | |||
| 86 |
1/2✓ Branch 1 taken 43724157 times.
✗ Branch 2 not taken.
|
78449750 | Helper::solveLocalSystem(this->fvGeometry(), handle, iv); |
| 87 | |||
| 88 | // make sure the inverse of A now carries zeroes in the zero rows, as | ||
| 89 | // well as CA and T. AB will have the zeroes already because we set the rows | ||
| 90 | // of B to zero above. | ||
| 91 |
4/4✓ Branch 0 taken 75879569 times.
✓ Branch 1 taken 43724157 times.
✓ Branch 2 taken 75879569 times.
✓ Branch 3 taken 43724157 times.
|
499599326 | for (const auto& zeroRowIndices : zeroRows) |
| 92 | { | ||
| 93 | 132125038 | const auto faceIdx = zeroRowIndices.second; | |
| 94 | 264250076 | A[zeroRowIndices.first] = 0.0; | |
| 95 | 396375114 | handle.CA()[faceIdx] = 0.0; | |
| 96 | 528500152 | handle.T()[faceIdx] = 0.0; | |
| 97 | |||
| 98 | // reset outside transmissibilities on surface grids | ||
| 99 | static constexpr int dim = IV::Traits::GridView::dimension; | ||
| 100 | static constexpr int dimWorld = IV::Traits::GridView::dimensionworld; | ||
| 101 | if constexpr (dim < dimWorld) | ||
| 102 | ✗ | std::for_each( handle.tijOutside()[faceIdx].begin(), | |
| 103 | ✗ | handle.tijOutside()[faceIdx].end(), | |
| 104 | ✗ | [] (auto& outsideTij) { outsideTij = 0.0; } ); | |
| 105 | } | ||
| 106 | } | ||
| 107 | 78449942 | } | |
| 108 | |||
| 109 | /*! | ||
| 110 | * \brief Assembles the vector of primary (cell) unknowns and (maybe) | ||
| 111 | * Dirichlet boundary conditions within an interaction volume. | ||
| 112 | * | ||
| 113 | * \param handle The data handle in which the vector is stored | ||
| 114 | * \param iv The mpfa-o interaction volume | ||
| 115 | * \param getU Lambda to obtain the desired cell/Dirichlet value from vol vars | ||
| 116 | */ | ||
| 117 | template< class DataHandle, class IV, class GetU > | ||
| 118 | 156906442 | void assembleU(DataHandle& handle, const IV& iv, const GetU& getU) | |
| 119 | { | ||
| 120 | 157140368 | auto& u = handle.uj(); | |
| 121 | 157140368 | Helper::resizeVector(u, iv.numKnowns()); | |
| 122 | |||
| 123 | // put the cell unknowns first, then Dirichlet values | ||
| 124 | 157140368 | typename IV::Traits::IndexSet::LocalIndexType i = 0; | |
| 125 |
4/4✓ Branch 0 taken 419417644 times.
✓ Branch 1 taken 95851318 times.
✓ Branch 2 taken 415647356 times.
✓ Branch 3 taken 94908746 times.
|
1620122782 | for (; i < iv.numScvs(); i++) |
| 126 | 1314098136 | u[i] = getU( this->elemVolVars()[iv.localScv(i).gridScvIndex()] ); | |
| 127 |
4/4✓ Branch 0 taken 2006852 times.
✓ Branch 1 taken 94908746 times.
✓ Branch 2 taken 2006852 times.
✓ Branch 3 taken 94908746 times.
|
634033984 | for (const auto& data : iv.dirichletData()) |
| 128 | 2736256 | u[i++] = getU( this->elemVolVars()[data.volVarIndex()] ); | |
| 129 | 156906442 | } | |
| 130 | |||
| 131 | private: | ||
| 132 | /*! | ||
| 133 | * \brief Assemble the matrices involved in the flux expressions | ||
| 134 | * across the scvfs inside an interaction volume as well as those involved | ||
| 135 | * in the interaction volume-local system of equations resulting from flux | ||
| 136 | * and solution continuity across the scvfs. | ||
| 137 | * | ||
| 138 | * Flux expressions: \f$\mathbf{f} = \mathbf{C} \bar{\mathbf{u}} + \mathbf{D} \mathbf{u}\f$. | ||
| 139 | * | ||
| 140 | * Continuity equations: \f$\mathbf{A} \, \bar{\mathbf{u}} = \mathbf{B} \, \mathbf{u}\f$. | ||
| 141 | * | ||
| 142 | * \note The matrices are expected to have been resized beforehand. | ||
| 143 | * | ||
| 144 | * \tparam IV The interaction volume type implementation | ||
| 145 | * \tparam TensorFunc Lambda to obtain the tensor w.r.t. which the local system is to be solved | ||
| 146 | * | ||
| 147 | * \param A The A matrix of the iv-local equation system | ||
| 148 | * \param B The B matrix of the iv-local equation system | ||
| 149 | * \param C The C matrix of the iv-local flux expressions | ||
| 150 | * \param D The D matrix of the iv-local flux expressions | ||
| 151 | * \param iv The mpfa-o interaction volume | ||
| 152 | * \param getT Lambda to evaluate the scv-wise tensors | ||
| 153 | * \param wijZeroThresh Threshold below which transmissibilities are taken to be zero. | ||
| 154 | * On the basis of this threshold, trivial (0 = 0) rows in the | ||
| 155 | * A matrix are identified and modified accordingly in order to | ||
| 156 | * avoid ending up with singular matrices. This can occur when the | ||
| 157 | * tensor is zero in some cells. | ||
| 158 | * \returns std::vector of pairs storing the local dof and the local face indices of those | ||
| 159 | * non-Dirichlet faces, for which the flux continuity condition results in a | ||
| 160 | * trivial 0 = 0 equation which could happen for zero tensors. | ||
| 161 | */ | ||
| 162 | template< class IV, class TensorFunc > | ||
| 163 | 43724349 | auto assembleLocalMatrices_(typename IV::Traits::MatVecTraits::AMatrix& A, | |
| 164 | typename IV::Traits::MatVecTraits::BMatrix& B, | ||
| 165 | typename IV::Traits::MatVecTraits::CMatrix& C, | ||
| 166 | typename IV::Traits::MatVecTraits::DMatrix& D, | ||
| 167 | typename IV::Traits::MatVecTraits::OmegaStorage& wijk, | ||
| 168 | IV& iv, const TensorFunc& getT, | ||
| 169 | Scalar<IV> wijZeroThresh) | ||
| 170 | { | ||
| 171 | using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; | ||
| 172 | static constexpr int dim = IV::Traits::GridView::dimension; | ||
| 173 | static constexpr int dimWorld = IV::Traits::GridView::dimensionworld; | ||
| 174 | |||
| 175 |
1/2✓ Branch 1 taken 8998756 times.
✗ Branch 2 not taken.
|
43724349 | std::vector< std::pair<LocalIndexType, LocalIndexType> > faceMarkers; |
| 176 |
1/2✓ Branch 1 taken 8998756 times.
✗ Branch 2 not taken.
|
43724349 | Helper::resizeVector(wijk, iv.numFaces()); |
| 177 | |||
| 178 | // if only Dirichlet faces are present in the iv, | ||
| 179 | // the matrices A, B & C are undefined and D = T | ||
| 180 |
2/2✓ Branch 0 taken 192 times.
✓ Branch 1 taken 8998564 times.
|
43017420 | if (iv.numUnknowns() == 0) |
| 181 | { | ||
| 182 | // resize & reset D matrix | ||
| 183 |
1/2✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
|
192 | Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0; |
| 184 | |||
| 185 | // Loop over all the faces, in this case these are all dirichlet boundaries | ||
| 186 |
2/2✓ Branch 0 taken 384 times.
✓ Branch 1 taken 192 times.
|
576 | for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) |
| 187 | { | ||
| 188 | 384 | const auto& curLocalScvf = iv.localScvf(faceIdx); | |
| 189 | 384 | const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); | |
| 190 | 384 | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); | |
| 191 | |||
| 192 | // get tensor in "positive" sub volume | ||
| 193 |
0/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
768 | const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); |
| 194 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
384 | const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex()); |
| 195 | 384 | const auto& posVolVars = this->elemVolVars()[posGlobalScv]; | |
| 196 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
384 | const auto tensor = getT(posVolVars); |
| 197 | |||
| 198 | // the omega factors of the "positive" sub volume | ||
| 199 |
0/4✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
768 | Helper::resizeVector(wijk[faceIdx], /*no outside scvs present*/1); |
| 200 | 384 | wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); | |
| 201 | |||
| 202 | 384 | const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); | |
| 203 |
2/2✓ Branch 0 taken 768 times.
✓ Branch 1 taken 384 times.
|
1152 | for (LocalIndexType localDir = 0; localDir < dim; localDir++) |
| 204 | { | ||
| 205 | 768 | const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) ); | |
| 206 | 768 | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); | |
| 207 | 3072 | D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; | |
| 208 | 3840 | D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir]; | |
| 209 | } | ||
| 210 | } | ||
| 211 | } | ||
| 212 | else | ||
| 213 | { | ||
| 214 | // resize & reset matrices | ||
| 215 |
1/2✓ Branch 1 taken 8998564 times.
✗ Branch 2 not taken.
|
43724157 | Helper::resizeMatrix(A, iv.numUnknowns(), iv.numUnknowns()); A = 0.0; |
| 216 |
1/2✓ Branch 1 taken 8998564 times.
✗ Branch 2 not taken.
|
43017228 | Helper::resizeMatrix(B, iv.numUnknowns(), iv.numKnowns()); B = 0.0; |
| 217 |
1/2✓ Branch 1 taken 8998564 times.
✗ Branch 2 not taken.
|
43017228 | Helper::resizeMatrix(C, iv.numFaces(), iv.numUnknowns()); C = 0.0; |
| 218 |
1/2✓ Branch 1 taken 8998564 times.
✗ Branch 2 not taken.
|
43017228 | Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0; |
| 219 | |||
| 220 |
2/2✓ Branch 0 taken 34527236 times.
✓ Branch 1 taken 8998564 times.
|
214330369 | for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) |
| 221 | { | ||
| 222 | 170606212 | const auto& curLocalScvf = iv.localScvf(faceIdx); | |
| 223 | 170606212 | const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); | |
| 224 | 170606212 | const auto curIsDirichlet = curLocalScvf.isDirichlet(); | |
| 225 | 170606212 | const auto curLocalDofIdx = curLocalScvf.localDofIndex(); | |
| 226 | |||
| 227 | // get tensor in "positive" sub volume | ||
| 228 | 170606212 | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); | |
| 229 |
4/4✓ Branch 0 taken 1775990 times.
✓ Branch 1 taken 5078860 times.
✓ Branch 2 taken 1775990 times.
✓ Branch 3 taken 5078860 times.
|
339068718 | const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); |
| 230 |
2/2✓ Branch 0 taken 1775990 times.
✓ Branch 1 taken 5078860 times.
|
170606212 | const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex()); |
| 231 | 170606212 | const auto& posVolVars = this->elemVolVars()[posGlobalScv]; | |
| 232 |
1/2✓ Branch 1 taken 2143706 times.
✗ Branch 2 not taken.
|
170606212 | const auto tensor = getT(posVolVars); |
| 233 | |||
| 234 | // the omega factors of the "positive" sub volume | ||
| 235 |
3/6✓ Branch 1 taken 2143706 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2143706 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2143706 times.
✗ Branch 8 not taken.
|
343356130 | Helper::resizeVector(wijk[faceIdx], neighborScvIndices.size()); |
| 236 | 170606212 | wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); | |
| 237 | |||
| 238 | using std::abs; | ||
| 239 | 170606212 | bool insideZeroWij = false; | |
| 240 | |||
| 241 | // go over the coordinate directions in the positive sub volume | ||
| 242 |
2/2✓ Branch 0 taken 69054472 times.
✓ Branch 1 taken 34527236 times.
|
511818636 | for (unsigned int localDir = 0; localDir < dim; localDir++) |
| 243 | { | ||
| 244 |
2/2✓ Branch 1 taken 34527236 times.
✓ Branch 2 taken 34527236 times.
|
341212424 | const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) ); |
| 245 | 341212424 | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); | |
| 246 | |||
| 247 |
2/2✓ Branch 0 taken 34527236 times.
✓ Branch 1 taken 34527236 times.
|
341212424 | if (otherLocalDofIdx == curLocalDofIdx) |
| 248 | { | ||
| 249 |
10/10✓ Branch 0 taken 19582696 times.
✓ Branch 1 taken 14944540 times.
✓ Branch 2 taken 19582696 times.
✓ Branch 3 taken 14944540 times.
✓ Branch 4 taken 19582696 times.
✓ Branch 5 taken 14944540 times.
✓ Branch 6 taken 19582696 times.
✓ Branch 7 taken 14944540 times.
✓ Branch 8 taken 19582696 times.
✓ Branch 9 taken 12800834 times.
|
850887354 | if (abs(wijk[faceIdx][0][localDir]) <= wijZeroThresh) |
| 250 | { | ||
| 251 |
2/2✓ Branch 0 taken 19421816 times.
✓ Branch 1 taken 160880 times.
|
75774286 | if (!curIsDirichlet) |
| 252 | { | ||
| 253 | 75436842 | insideZeroWij = true; | |
| 254 |
2/4✓ Branch 1 taken 19421816 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19421816 times.
✗ Branch 5 not taken.
|
150873684 | faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) ); |
| 255 | } | ||
| 256 | } | ||
| 257 | } | ||
| 258 | |||
| 259 | // if we are not on a Dirichlet face, add entries associated with unknown face pressures | ||
| 260 | // i.e. in matrix C and maybe A (if current face is not a Dirichlet face) | ||
| 261 |
2/2✓ Branch 0 taken 68257046 times.
✓ Branch 1 taken 797426 times.
|
341212424 | if (!otherLocalScvf.isDirichlet()) |
| 262 | { | ||
| 263 |
10/10✓ Branch 0 taken 67716400 times.
✓ Branch 1 taken 540646 times.
✓ Branch 2 taken 67716400 times.
✓ Branch 3 taken 540646 times.
✓ Branch 4 taken 67716400 times.
✓ Branch 5 taken 540646 times.
✓ Branch 6 taken 67716400 times.
✓ Branch 7 taken 540646 times.
✓ Branch 8 taken 63765996 times.
✓ Branch 9 taken 400150 times.
|
1694614010 | C[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; |
| 264 |
2/2✓ Branch 0 taken 67716400 times.
✓ Branch 1 taken 540646 times.
|
339740982 | if (!curIsDirichlet) |
| 265 | 1695495508 | A[curLocalDofIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; | |
| 266 | } | ||
| 267 | // the current face is a Dirichlet face and creates entries in D & maybe B | ||
| 268 | else | ||
| 269 | { | ||
| 270 |
10/10✓ Branch 0 taken 256780 times.
✓ Branch 1 taken 540646 times.
✓ Branch 2 taken 256780 times.
✓ Branch 3 taken 540646 times.
✓ Branch 4 taken 256780 times.
✓ Branch 5 taken 540646 times.
✓ Branch 6 taken 256780 times.
✓ Branch 7 taken 540646 times.
✓ Branch 8 taken 200764 times.
✓ Branch 9 taken 400150 times.
|
7160698 | D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; |
| 271 |
2/2✓ Branch 0 taken 256780 times.
✓ Branch 1 taken 540646 times.
|
1471442 | if (!curIsDirichlet) |
| 272 | 2386764 | B[curLocalDofIdx][otherLocalDofIdx] += wijk[faceIdx][0][localDir]; | |
| 273 | } | ||
| 274 | |||
| 275 | // add entries related to pressures at the scv centers (dofs) | ||
| 276 | 341212424 | const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); | |
| 277 |
10/10✓ Branch 0 taken 67973180 times.
✓ Branch 1 taken 1081292 times.
✓ Branch 2 taken 67973180 times.
✓ Branch 3 taken 1081292 times.
✓ Branch 4 taken 67973180 times.
✓ Branch 5 taken 1081292 times.
✓ Branch 6 taken 67973180 times.
✓ Branch 7 taken 1081292 times.
✓ Branch 8 taken 63966760 times.
✓ Branch 9 taken 800300 times.
|
1701774708 | D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir]; |
| 278 | |||
| 279 |
2/2✓ Branch 0 taken 67973180 times.
✓ Branch 1 taken 1081292 times.
|
341212424 | if (!curIsDirichlet) |
| 280 | 1697882272 | B[curLocalDofIdx][posScvLocalDofIdx] -= wijk[faceIdx][0][localDir]; | |
| 281 | } | ||
| 282 | |||
| 283 | // If we are on an interior face, add values from negative sub volume | ||
| 284 |
2/2✓ Branch 0 taken 30343664 times.
✓ Branch 1 taken 4183572 times.
|
170606212 | if (!curGlobalScvf.boundary()) |
| 285 | { | ||
| 286 | // loop over all the outside neighbors of this face and add entries | ||
| 287 |
4/4✓ Branch 0 taken 30384580 times.
✓ Branch 1 taken 30343664 times.
✓ Branch 2 taken 30384580 times.
✓ Branch 3 taken 30343664 times.
|
476327588 | for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) |
| 288 | { | ||
| 289 | 158803140 | const auto idxOnScvf = idxInOutside+1; | |
| 290 |
4/4✓ Branch 0 taken 1705182 times.
✓ Branch 1 taken 5006592 times.
✓ Branch 2 taken 1705182 times.
✓ Branch 3 taken 5006592 times.
|
317606280 | const auto& negLocalScv = iv.localScv( neighborScvIndices[idxOnScvf] ); |
| 291 |
2/2✓ Branch 0 taken 1705182 times.
✓ Branch 1 taken 5006592 times.
|
158803140 | const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.gridScvIndex()); |
| 292 | 158803140 | const auto& negVolVars = this->elemVolVars()[negGlobalScv]; | |
| 293 | 158803140 | const auto negTensor = getT(negVolVars); | |
| 294 | |||
| 295 | // On surface grids, use outside face for "negative" transmissibility calculation | ||
| 296 | 158803140 | const auto& scvf = dim < dimWorld ? this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside) | |
| 297 | : curGlobalScvf; | ||
| 298 | 158803140 | wijk[faceIdx][idxOnScvf] = computeMpfaTransmissibility<EG>(this->fvGeometry(), negLocalScv, scvf, negTensor, negVolVars.extrusionFactor()); | |
| 299 | |||
| 300 | // flip sign on surface grids (since we used the "outside" normal) | ||
| 301 | if (dim < dimWorld) | ||
| 302 | 1354702 | wijk[faceIdx][idxOnScvf] *= -1.0; | |
| 303 | |||
| 304 | // go over the coordinate directions in the negative sub volume | ||
| 305 |
2/2✓ Branch 0 taken 60769160 times.
✓ Branch 1 taken 30384580 times.
|
476409420 | for (int localDir = 0; localDir < dim; localDir++) |
| 306 | { | ||
| 307 | 317606280 | const auto otherLocalScvfIdx = negLocalScv.localScvfIndex(localDir); | |
| 308 |
2/2✓ Branch 0 taken 30384580 times.
✓ Branch 1 taken 30384580 times.
|
317606280 | const auto& otherLocalScvf = iv.localScvf(otherLocalScvfIdx); |
| 309 | 317606280 | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); | |
| 310 | |||
| 311 | // check for zero transmissibilities (skip if inside has been zero already) | ||
| 312 |
4/4✓ Branch 0 taken 30384580 times.
✓ Branch 1 taken 30384580 times.
✓ Branch 2 taken 13665548 times.
✓ Branch 3 taken 16719032 times.
|
317606280 | if (otherLocalDofIdx == curLocalDofIdx && !insideZeroWij) |
| 313 |
10/10✓ Branch 0 taken 212284 times.
✓ Branch 1 taken 13453264 times.
✓ Branch 2 taken 212284 times.
✓ Branch 3 taken 13453264 times.
✓ Branch 4 taken 212284 times.
✓ Branch 5 taken 13453264 times.
✓ Branch 6 taken 212284 times.
✓ Branch 7 taken 13453264 times.
✓ Branch 8 taken 212284 times.
✓ Branch 9 taken 13453264 times.
|
434181200 | if (abs(wijk[faceIdx][idxOnScvf][localDir]) <= wijZeroThresh) |
| 314 |
2/4✓ Branch 1 taken 212284 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 212284 times.
✗ Branch 5 not taken.
|
885454 | faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) ); |
| 315 | |||
| 316 |
2/2✓ Branch 0 taken 60485294 times.
✓ Branch 1 taken 283866 times.
|
317606280 | if (!otherLocalScvf.isDirichlet()) |
| 317 | 1902671700 | A[curLocalDofIdx][otherLocalDofIdx] += wijk[faceIdx][idxOnScvf][localDir]; | |
| 318 | else | ||
| 319 | 2965980 | B[curLocalDofIdx][otherLocalDofIdx] -= wijk[faceIdx][idxOnScvf][localDir]; | |
| 320 | |||
| 321 | // add entries to matrix B | ||
| 322 | 1905637680 | B[curLocalDofIdx][negLocalScv.localDofIndex()] += wijk[faceIdx][idxOnScvf][localDir]; | |
| 323 | } | ||
| 324 | } | ||
| 325 | } | ||
| 326 | } | ||
| 327 | } | ||
| 328 | |||
| 329 | 43724349 | return faceMarkers; | |
| 330 | } | ||
| 331 | }; | ||
| 332 | |||
| 333 | } // end namespace | ||
| 334 | |||
| 335 | #endif | ||
| 336 |