| 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-FileCopyrightText: 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 |
5/5✓ Branch 0 taken 398494 times.
✓ Branch 1 taken 47445360 times.
✓ Branch 3 taken 12065201 times.
✓ Branch 4 taken 53662 times.
✓ Branch 2 taken 3231066 times.
|
67649652 | 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 | 78926102 | void assembleMatrices(DataHandle& handle, IV& iv, const TensorFunc& getT, Scalar<IV> wijZeroThresh = 0.0) | |
| 63 | { | ||
| 64 | 78926102 | 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 43255308 times.
✓ Branch 1 taken 192 times.
|
77512244 | if (iv.numUnknowns() > 0) |
| 68 | { | ||
| 69 | // for now, B is carried in what will be AB after solve | ||
| 70 | 78925910 | auto& B = handle.AB(); | |
| 71 | 78925910 | 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 |
2/2✓ Branch 0 taken 75879569 times.
✓ Branch 1 taken 43962237 times.
|
211050948 | 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.
|
660810816 | for (auto& row : A) |
| 80 | 528685778 | row[zeroRowDofIdx] = 0.0; | |
| 81 | 132125038 | A[zeroRowDofIdx] = 0.0; | |
| 82 | 132125038 | A[zeroRowDofIdx][zeroRowDofIdx] = 1.0; | |
| 83 | 264250076 | B[zeroRowDofIdx] = 0.0; | |
| 84 | } | ||
| 85 | |||
| 86 |
1/2✓ Branch 1 taken 43962237 times.
✗ Branch 2 not taken.
|
78925910 | 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 |
2/2✓ Branch 0 taken 75879569 times.
✓ Branch 1 taken 43962237 times.
|
211050948 | for (const auto& zeroRowIndices : zeroRows) |
| 92 | { | ||
| 93 | 132125038 | const auto faceIdx = zeroRowIndices.second; | |
| 94 | 132125038 | A[zeroRowIndices.first] = 0.0; | |
| 95 | 132125038 | handle.CA()[faceIdx] = 0.0; | |
| 96 | 264250076 | 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 | 78926102 | } | |
| 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 | 157850454 | void assembleU(DataHandle& handle, const IV& iv, const GetU& getU) | |
| 119 | { | ||
| 120 | 157850454 | auto& u = handle.uj(); | |
| 121 | 155965310 | Helper::resizeVector(u, iv.numKnowns()); | |
| 122 | |||
| 123 | // put the cell unknowns first, then Dirichlet values | ||
| 124 | 157850454 | typename IV::Traits::IndexSet::LocalIndexType i = 0; | |
| 125 |
2/2✓ Branch 0 taken 419890828 times.
✓ Branch 1 taken 96089398 times.
|
816781594 | for (; i < iv.numScvs(); i++) |
| 126 | 658931140 | u[i] = getU( this->elemVolVars()[iv.localScv(i).gridScvIndex()] ); | |
| 127 |
2/2✓ Branch 0 taken 2009828 times.
✓ Branch 1 taken 95146826 times.
|
160592662 | for (const auto& data : iv.dirichletData()) |
| 128 | 2742208 | u[i++] = getU( this->elemVolVars()[data.volVarIndex()] ); | |
| 129 | 157850454 | } | |
| 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 | 78926102 | 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 43255500 times.
✗ Branch 2 not taken.
|
78926102 | std::vector< std::pair<LocalIndexType, LocalIndexType> > faceMarkers; |
| 176 |
3/4✓ Branch 1 taken 43255500 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 192 times.
✓ Branch 4 taken 43255308 times.
|
155024488 | 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 43255308 times.
|
77512244 | if (iv.numUnknowns() == 0) |
| 181 | { | ||
| 182 | // resize & reset D matrix | ||
| 183 |
1/2✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
|
384 | 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 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
384 | const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); |
| 190 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
384 | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); |
| 191 | |||
| 192 | // get tensor in "positive" sub volume | ||
| 193 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
384 | 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 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
384 | const auto& posVolVars = this->elemVolVars()[posGlobalScv]; |
| 196 | 384 | const auto tensor = getT(posVolVars); | |
| 197 | |||
| 198 | // the omega factors of the "positive" sub volume | ||
| 199 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
384 | 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 | 768 | D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; | |
| 208 | 768 | D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir]; | |
| 209 | } | ||
| 210 | } | ||
| 211 | } | ||
| 212 | else | ||
| 213 | { | ||
| 214 | // resize & reset matrices | ||
| 215 |
3/5✓ Branch 0 taken 2827716 times.
✓ Branch 1 taken 43962237 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43255308 times.
✗ Branch 5 not taken.
|
241019304 | Helper::resizeMatrix(A, iv.numUnknowns(), iv.numUnknowns()); A = 0.0; |
| 216 |
3/5✓ Branch 0 taken 2827716 times.
✓ Branch 1 taken 43962237 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43255308 times.
✗ Branch 5 not taken.
|
239605446 | Helper::resizeMatrix(B, iv.numUnknowns(), iv.numKnowns()); B = 0.0; |
| 217 |
2/3✓ Branch 0 taken 2827716 times.
✓ Branch 1 taken 43962237 times.
✗ Branch 2 not taken.
|
162093394 | Helper::resizeMatrix(C, iv.numFaces(), iv.numUnknowns()); C = 0.0; |
| 218 |
1/2✓ Branch 1 taken 43255308 times.
✗ Branch 2 not taken.
|
162093394 | Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0; |
| 219 | |||
| 220 |
2/2✓ Branch 0 taken 171317476 times.
✓ Branch 1 taken 43962237 times.
|
387033626 | for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) |
| 221 | { | ||
| 222 | 308107716 | const auto& curLocalScvf = iv.localScvf(faceIdx); | |
| 223 |
2/2✓ Branch 1 taken 5385086 times.
✓ Branch 2 taken 9383860 times.
|
308107716 | const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); |
| 224 |
2/2✓ Branch 0 taken 5385086 times.
✓ Branch 1 taken 9383860 times.
|
308107716 | const auto curIsDirichlet = curLocalScvf.isDirichlet(); |
| 225 |
2/2✓ Branch 0 taken 5385086 times.
✓ Branch 1 taken 9383860 times.
|
308107716 | const auto curLocalDofIdx = curLocalScvf.localDofIndex(); |
| 226 | |||
| 227 | // get tensor in "positive" sub volume | ||
| 228 |
2/2✓ Branch 0 taken 5385086 times.
✓ Branch 1 taken 9383860 times.
|
308107716 | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); |
| 229 |
2/2✓ Branch 0 taken 5385086 times.
✓ Branch 1 taken 9383860 times.
|
308107716 | const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); |
| 230 |
2/2✓ Branch 0 taken 5385086 times.
✓ Branch 1 taken 9383860 times.
|
308107716 | const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex()); |
| 231 |
1/2✓ Branch 2 taken 2143706 times.
✗ Branch 3 not taken.
|
308107716 | const auto& posVolVars = this->elemVolVars()[posGlobalScv]; |
| 232 |
1/2✓ Branch 1 taken 2143706 times.
✗ Branch 2 not taken.
|
308107716 | const auto tensor = getT(posVolVars); |
| 233 | |||
| 234 | // the omega factors of the "positive" sub volume | ||
| 235 |
1/2✓ Branch 1 taken 2143706 times.
✗ Branch 2 not taken.
|
308107716 | Helper::resizeVector(wijk[faceIdx], neighborScvIndices.size()); |
| 236 | 308107716 | wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); | |
| 237 | |||
| 238 | using std::abs; | ||
| 239 | 308107716 | bool insideZeroWij = false; | |
| 240 | |||
| 241 | // go over the coordinate directions in the positive sub volume | ||
| 242 |
2/2✓ Branch 0 taken 342634952 times.
✓ Branch 1 taken 171317476 times.
|
924323148 | for (unsigned int localDir = 0; localDir < dim; localDir++) |
| 243 | { | ||
| 244 |
2/2✓ Branch 1 taken 171317476 times.
✓ Branch 2 taken 171317476 times.
|
616215432 | const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) ); |
| 245 | 616215432 | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); | |
| 246 | |||
| 247 |
2/2✓ Branch 0 taken 171317476 times.
✓ Branch 1 taken 171317476 times.
|
616215432 | if (otherLocalDofIdx == curLocalDofIdx) |
| 248 | { | ||
| 249 |
2/2✓ Branch 0 taken 75774286 times.
✓ Branch 1 taken 95543190 times.
|
308107716 | if (abs(wijk[faceIdx][0][localDir]) <= wijZeroThresh) |
| 250 | { | ||
| 251 |
2/2✓ Branch 0 taken 75436842 times.
✓ Branch 1 taken 337444 times.
|
131965876 | if (!curIsDirichlet) |
| 252 | { | ||
| 253 |
1/2✓ Branch 1 taken 75436842 times.
✗ Branch 2 not taken.
|
131451868 | insideZeroWij = true; |
| 254 |
1/2✓ Branch 1 taken 75436842 times.
✗ Branch 2 not taken.
|
131451868 | 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 341157558 times.
✓ Branch 1 taken 1477394 times.
|
616215432 | if (!otherLocalScvf.isDirichlet()) |
| 262 | { | ||
| 263 |
2/2✓ Branch 0 taken 340171696 times.
✓ Branch 1 taken 985862 times.
|
614058070 | C[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; |
| 264 |
2/2✓ Branch 0 taken 340171696 times.
✓ Branch 1 taken 985862 times.
|
614058070 | if (!curIsDirichlet) |
| 265 | 612626992 | 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 |
2/2✓ Branch 0 taken 491532 times.
✓ Branch 1 taken 985862 times.
|
2157362 | D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; |
| 271 |
2/2✓ Branch 0 taken 491532 times.
✓ Branch 1 taken 985862 times.
|
2157362 | if (!curIsDirichlet) |
| 272 | 726284 | B[curLocalDofIdx][otherLocalDofIdx] += wijk[faceIdx][0][localDir]; | |
| 273 | } | ||
| 274 | |||
| 275 | // add entries related to pressures at the scv centers (dofs) | ||
| 276 | 616215432 | const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); | |
| 277 |
2/2✓ Branch 0 taken 340663228 times.
✓ Branch 1 taken 1971724 times.
|
616215432 | D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir]; |
| 278 | |||
| 279 |
2/2✓ Branch 0 taken 340663228 times.
✓ Branch 1 taken 1971724 times.
|
616215432 | if (!curIsDirichlet) |
| 280 | 613353276 | 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 158997328 times.
✓ Branch 1 taken 12320148 times.
|
308107716 | if (!curGlobalScvf.boundary()) |
| 285 | { | ||
| 286 | // loop over all the outside neighbors of this face and add entries | ||
| 287 |
2/2✓ Branch 0 taken 159038244 times.
✓ Branch 1 taken 158997328 times.
|
575342900 | for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) |
| 288 | { | ||
| 289 | 287691908 | const auto idxOnScvf = idxInOutside+1; | |
| 290 |
2/2✓ Branch 0 taken 3140182 times.
✓ Branch 1 taken 7178600 times.
|
287691908 | const auto& negLocalScv = iv.localScv( neighborScvIndices[idxOnScvf] ); |
| 291 |
2/2✓ Branch 0 taken 3140182 times.
✓ Branch 1 taken 7178600 times.
|
287691908 | const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.gridScvIndex()); |
| 292 | 287691908 | const auto& negVolVars = this->elemVolVars()[negGlobalScv]; | |
| 293 | 287691908 | const auto negTensor = getT(negVolVars); | |
| 294 | |||
| 295 | // On surface grids, use outside face for "negative" transmissibility calculation | ||
| 296 | 287691908 | const auto& scvf = dim < dimWorld ? this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside) | |
| 297 | : curGlobalScvf; | ||
| 298 | 287691908 | 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 318076488 times.
✓ Branch 1 taken 159038244 times.
|
863075724 | for (int localDir = 0; localDir < dim; localDir++) |
| 306 | { | ||
| 307 |
2/2✓ Branch 1 taken 159038244 times.
✓ Branch 2 taken 159038244 times.
|
575383816 | const auto otherLocalScvfIdx = negLocalScv.localScvfIndex(localDir); |
| 308 |
2/2✓ Branch 0 taken 159038244 times.
✓ Branch 1 taken 159038244 times.
|
575383816 | const auto& otherLocalScvf = iv.localScvf(otherLocalScvfIdx); |
| 309 | 575383816 | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); | |
| 310 | |||
| 311 | // check for zero transmissibilities (skip if inside has been zero already) | ||
| 312 |
4/4✓ Branch 0 taken 159038244 times.
✓ Branch 1 taken 159038244 times.
✓ Branch 2 taken 87071344 times.
✓ Branch 3 taken 71966900 times.
|
575383816 | if (otherLocalDofIdx == curLocalDofIdx && !insideZeroWij) |
| 313 |
2/2✓ Branch 0 taken 442727 times.
✓ Branch 1 taken 86628617 times.
|
160477140 | if (abs(wijk[faceIdx][idxOnScvf][localDir]) <= wijZeroThresh) |
| 314 |
1/2✓ Branch 1 taken 442727 times.
✗ Branch 2 not taken.
|
673170 | faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) ); |
| 315 | |||
| 316 |
2/2✓ Branch 0 taken 317582158 times.
✓ Branch 1 taken 494330 times.
|
575383816 | if (!otherLocalScvf.isDirichlet()) |
| 317 | 574679022 | A[curLocalDofIdx][otherLocalDofIdx] += wijk[faceIdx][idxOnScvf][localDir]; | |
| 318 | else | ||
| 319 | 704794 | B[curLocalDofIdx][otherLocalDofIdx] -= wijk[faceIdx][idxOnScvf][localDir]; | |
| 320 | |||
| 321 | // add entries to matrix B | ||
| 322 | 575383816 | B[curLocalDofIdx][negLocalScv.localDofIndex()] += wijk[faceIdx][idxOnScvf][localDir]; | |
| 323 | } | ||
| 324 | } | ||
| 325 | } | ||
| 326 | } | ||
| 327 | } | ||
| 328 | |||
| 329 | 78926102 | return faceMarkers; | |
| 330 | ✗ | } | |
| 331 | }; | ||
| 332 | |||
| 333 | } // end namespace | ||
| 334 | |||
| 335 | #endif | ||
| 336 |