diff --git a/CHANGELOG.md b/CHANGELOG.md index bd94e7c656805179328a2f51476940b3a04aa988..2499c7d10a66f10ff0368123b0eacfd60cef5407 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -41,6 +41,9 @@ Then code possibly fails to compile. The fix is to implement the same caching co - __Box__: `scv/scvf.geometry()` have been deprecated, use `fvGeometry.geometry(scv/scvf)` instead - __Box__: `scv/scvf.corner(i)` have been deprecated, use `fvGeometry.geometry(scv/scvf).corner(i)` instead +- __Cell centered__: The `computeTpfa/MpfaTransmissibilities()` interfaces now require an additional parameter `fvGeometry` +- __Staggered__: `fluxVars.advectivFluxForCellCenter()/inflowOutflowBoundaryFlux()` interfaces now require parameter `fvGeometry` instead of `element` +- __Discretization__: `Extrusion::area/volume(scvf/scv)` have been deprecated, use `Extrusion::area/volume(fvGeometry, scvf/scv)` instead ### New experimental features (possibly subject to backwards-incompatible changes in the future) diff --git a/dumux/assembly/cclocalresidual.hh b/dumux/assembly/cclocalresidual.hh index e69b73f760260a2c45b1e3fade473e45c59c03ae..e02236c6b9aef48bf7f6a9a9e064a0fcf043fcf7 100644 --- a/dumux/assembly/cclocalresidual.hh +++ b/dumux/assembly/cclocalresidual.hh @@ -104,7 +104,7 @@ public: // multiply neumann fluxes with the area and the extrusion factor const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); - neumannFluxes *= Extrusion::area(scvf)*elemVolVars[scv].extrusionFactor(); + neumannFluxes *= Extrusion::area(fvGeometry, scvf)*elemVolVars[scv].extrusionFactor(); flux += neumannFluxes; } diff --git a/dumux/assembly/cvfelocalresidual.hh b/dumux/assembly/cvfelocalresidual.hh index 00122079f3c50fdc2c72466478e49c75401eea3f..a360aa32e51e0f601a49dfca12969ce82e4518c0 100644 --- a/dumux/assembly/cvfelocalresidual.hh +++ b/dumux/assembly/cvfelocalresidual.hh @@ -155,7 +155,7 @@ public: auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); // multiply neumann fluxes with the area and the extrusion factor - neumannFluxes *= Extrusion::area(scvf)*elemVolVars[scv].extrusionFactor(); + neumannFluxes *= Extrusion::area(fvGeometry, scvf)*elemVolVars[scv].extrusionFactor(); // only add fluxes to equations for which Neumann is set for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx) @@ -204,7 +204,7 @@ public: storage *= curVolVars.extrusionFactor(); storage -= prevStorage; - storage *= Extrusion::volume(scv); + storage *= Extrusion::volume(fvGeometry, scv); storage /= this->timeLoop().timeStepSize(); residual[scv.localDofIndex()] += storage; diff --git a/dumux/assembly/fclocalresidual.hh b/dumux/assembly/fclocalresidual.hh index 5d4ff9051f0b1c0ecef361758c9d0173f897f39a..15a2cca28e15afcbd95e25aabaed07b50dcb1f55 100644 --- a/dumux/assembly/fclocalresidual.hh +++ b/dumux/assembly/fclocalresidual.hh @@ -167,7 +167,7 @@ public: storage *= curVolVars.extrusionFactor(); storage -= prevStorage; - storage *= Extrusion::volume(scv); + storage *= Extrusion::volume(fvGeometry, scv); storage /= this->timeLoop().timeStepSize(); residual[scv.localDofIndex()] += storage; diff --git a/dumux/assembly/fvlocalresidual.hh b/dumux/assembly/fvlocalresidual.hh index 55872ff55e10aa9ca1e4ab30c56cb71c12143e2a..b5d0fea00498a229146285e27a6a7608a8f5f471 100644 --- a/dumux/assembly/fvlocalresidual.hh +++ b/dumux/assembly/fvlocalresidual.hh @@ -114,7 +114,7 @@ public: auto localScvIdx = scv.localDofIndex(); const auto& volVars = elemVolVars[scv]; storage[localScvIdx] = asImp().computeStorage(problem, scv, volVars); - storage[localScvIdx] *= Extrusion::volume(scv) * volVars.extrusionFactor(); + storage[localScvIdx] *= Extrusion::volume(fvGeometry, scv) * volVars.extrusionFactor(); } return storage; @@ -307,7 +307,7 @@ public: storage *= curVolVars.extrusionFactor(); storage -= prevStorage; - storage *= Extrusion::volume(scv); + storage *= Extrusion::volume(fvGeometry, scv); storage /= timeLoop_->timeStepSize(); residual[scv.localDofIndex()] += storage; @@ -336,7 +336,7 @@ public: //! Compute source with the model specific storage residual const auto& curVolVars = curElemVolVars[scv]; NumEqVector source = asImp().computeSource(problem, element, fvGeometry, curElemVolVars, scv); - source *= Extrusion::volume(scv)*curVolVars.extrusionFactor(); + source *= Extrusion::volume(fvGeometry, scv)*curVolVars.extrusionFactor(); //! subtract source from local rate (sign convention in user interface) residual[scv.localDofIndex()] -= source; diff --git a/dumux/assembly/staggeredlocalresidual.hh b/dumux/assembly/staggeredlocalresidual.hh index 927f01d14412c382d718aa5b7dd31f13ea7d5fff..15826b6bbf361e20c44593e7a2aad525cb6858f2 100644 --- a/dumux/assembly/staggeredlocalresidual.hh +++ b/dumux/assembly/staggeredlocalresidual.hh @@ -125,7 +125,7 @@ public: // subtract the source term from the local rate auto source = asImp_().computeSourceForCellCenter(problem, element, fvGeometry, curElemVolVars, curElemFaceVars, scv); - source *= Extrusion::volume(scv)*curExtrusionFactor; + source *= Extrusion::volume(fvGeometry, scv)*curExtrusionFactor; residual -= source; } @@ -171,7 +171,7 @@ public: storage = std::move(curCCStorage); storage -= std::move(prevCCStorage); - storage *= Extrusion::volume(scv); + storage *= Extrusion::volume(fvGeometry, scv); storage /= timeLoop_->timeStepSize(); residual += storage; @@ -250,7 +250,7 @@ public: faceScvCenter *= 0.5; FaceSubControlVolume faceScv(faceScvCenter, 0.5*scv.volume()); - source *= Extrusion::volume(faceScv)*extrusionFactor; + source *= Extrusion::volume(fvGeometry, faceScv)*extrusionFactor; residual -= source; } @@ -292,7 +292,7 @@ public: faceScvCenter *= 0.5; FaceSubControlVolume faceScv(faceScvCenter, 0.5*scv.volume()); - storage *= Extrusion::volume(faceScv)*extrusionFactor; + storage *= Extrusion::volume(fvGeometry, faceScv)*extrusionFactor; storage /= timeLoop_->timeStepSize(); residual += storage; diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index bf040be778207f63a132af5dd51a13f460cfadd4..a814d2a49c1c8b6c8a10a2cdf07fcbfb243aa080 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -444,7 +444,7 @@ public: // Add the contributions to the dof source values // We divide by the volume. In the local residual this will be multiplied with the same // factor again. That's because the user specifies absolute values in kg/s. - const auto volume = Extrusion::volume(scv)*elemVolVars[scv].extrusionFactor(); + const auto volume = Extrusion::volume(fvGeometry, scv)*elemVolVars[scv].extrusionFactor(); for (const auto& ps : pointSourceMap_.at(key)) { diff --git a/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh b/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh index 0159500ac746b95481aa46fb2c411fa8b53ced41..ed0ddb0a4401259130588e2d9644015375bcc4ce 100644 --- a/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh +++ b/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh @@ -47,6 +47,7 @@ namespace Dumux { * \param extrusionFactor The extrusion factor of the scv */ template +[[deprecated("Will be removed after release 3.6. Use interface with additional fvGeometry parameter instead.")]] Dune::FieldVector computeMpfaTransmissibility(const IVSubControlVolume& scv, const typename EG::SubControlVolumeFace& scvf, @@ -64,6 +65,38 @@ computeMpfaTransmissibility(const IVSubControlVolume& scv, return wijk; } +/*! + * \ingroup CCMpfaDiscretization + * \brief Free function to evaluate the Mpfa transmissibility associated + * with the flux (in the form of flux = t*gradU) across a + * sub-control volume face stemming from a given sub-control + * volume with corresponding tensor t. + * + * \param fvGeometry The element-centered control volume geometry + * \param scv The iv-local sub-control volume + * \param scvf The grid sub-control volume face + * \param t The tensor living in the scv + * \param extrusionFactor The extrusion factor of the scv + */ +template +Dune::FieldVector +computeMpfaTransmissibility(const EG& fvGeometry, + const IVSubControlVolume& scv, + const typename EG::SubControlVolumeFace& scvf, + const Tensor& t, + typename IVSubControlVolume::ctype extrusionFactor) +{ + Dune::FieldVector wijk; + for (unsigned int dir = 0; dir < IVSubControlVolume::myDimension; ++dir) + wijk[dir] = vtmv(scvf.unitOuterNormal(), t, scv.nu(dir)); + + using Extrusion = Extrusion_t; + wijk *= Extrusion::area(fvGeometry, scvf)*extrusionFactor; + wijk /= scv.detX(); + + return wijk; +} + /*! * \ingroup CCMpfaDiscretization * \brief Free function to evaluate the Mpfa transmissibility associated @@ -77,6 +110,7 @@ computeMpfaTransmissibility(const IVSubControlVolume& scv, * \param extrusionFactor The extrusion factor of the scv */ template< class EG, class IVSubControlVolume, class Tensor, std::enable_if_t< Dune::IsNumber::value, int > = 1 > +[[deprecated("Will be removed after release 3.6. Use interface with additional fvGeometry parameter instead.")]] Dune::FieldVector computeMpfaTransmissibility(const IVSubControlVolume& scv, const typename EG::SubControlVolumeFace& scvf, @@ -94,6 +128,38 @@ computeMpfaTransmissibility(const IVSubControlVolume& scv, return wijk; } +/*! + * \ingroup CCMpfaDiscretization + * \brief Free function to evaluate the Mpfa transmissibility associated + * with the flux (in the form of flux = t*gradU) across a + * sub-control volume face stemming from a given sub-control + * volume with corresponding tensor t, where t is a scalar. + * + * \param fvGeometry The element-centered control volume geometry + * \param scv The iv-local sub-control volume + * \param scvf The grid sub-control volume face + * \param t the scalar quantity living in the scv + * \param extrusionFactor The extrusion factor of the scv + */ +template< class EG, class IVSubControlVolume, class Tensor, std::enable_if_t< Dune::IsNumber::value, int > = 1 > +Dune::FieldVector +computeMpfaTransmissibility(const EG& fvGeometry, + const IVSubControlVolume& scv, + const typename EG::SubControlVolumeFace& scvf, + const Tensor& t, + typename IVSubControlVolume::ctype extrusionFactor) +{ + Dune::FieldVector wijk; + for (unsigned int dir = 0; dir < IVSubControlVolume::myDimension; ++dir) + wijk[dir] = vtmv(scvf.unitOuterNormal(), t, scv.nu(dir)); + + using Extrusion = Extrusion_t; + wijk *= Extrusion::area(fvGeometry, scvf)*extrusionFactor; + wijk /= scv.detX(); + + return wijk; +} + } // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh index af558765ae7e85b90c340469ebaa73ad85454f16..1e8d13f80af9207ac0ad91cba3ca600cc88b482b 100644 --- a/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh +++ b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh @@ -214,20 +214,20 @@ class InteractionVolumeAssemblerBase rho /= numOutsideFaces + 1; deltaG[localDofIdx] -= alpha_inside; - deltaG[localDofIdx] *= rho*Extrusion::area(curGlobalScvf); + deltaG[localDofIdx] *= rho*Extrusion::area(fvGeometry(), curGlobalScvf); } // use density resulting from Dirichlet BCs else rho = getRho(elemVolVars()[curGlobalScvf.outsideScvIdx()]); // add "inside" & "outside" alphas to gravity containers - g[faceIdx] = alpha_inside*rho*Extrusion::area(curGlobalScvf); + g[faceIdx] = alpha_inside*rho*Extrusion::area(fvGeometry(), curGlobalScvf); if (isSurfaceGrid) { unsigned int i = 0; for (const auto& alpha : alpha_outside) - outsideG[faceIdx][i++] = alpha*rho*Extrusion::area(curGlobalScvf); + outsideG[faceIdx][i++] = alpha*rho*Extrusion::area(fvGeometry(), curGlobalScvf); } } diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh index 5b0dc2ea74b7ea1e6ce4c64a4bf6511db234aeb6..0e954a09715786faaaded2a84b09221b7c948241 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh @@ -209,7 +209,7 @@ private: // the omega factors of the "positive" sub volume Helper::resizeVector(wijk[faceIdx], /*no outside scvs present*/1); - wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); + wijk[faceIdx][0] = computeMpfaTransmissibility(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); for (LocalIndexType localDir = 0; localDir < dim; localDir++) @@ -245,7 +245,7 @@ private: // the omega factors of the "positive" sub volume Helper::resizeVector(wijk[faceIdx], neighborScvIndices.size()); - wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); + wijk[faceIdx][0] = computeMpfaTransmissibility(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); using std::abs; bool insideZeroWij = false; @@ -307,7 +307,7 @@ private: // On surface grids, use outside face for "negative" transmissibility calculation const auto& scvf = dim < dimWorld ? this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside) : curGlobalScvf; - wijk[faceIdx][idxOnScvf] = computeMpfaTransmissibility(negLocalScv, scvf, negTensor, negVolVars.extrusionFactor()); + wijk[faceIdx][idxOnScvf] = computeMpfaTransmissibility(this->fvGeometry(), negLocalScv, scvf, negTensor, negVolVars.extrusionFactor()); // flip sign on surface grids (since we used the "outside" normal) if (dim < dimWorld) diff --git a/dumux/discretization/cellcentered/tpfa/computetransmissibility.hh b/dumux/discretization/cellcentered/tpfa/computetransmissibility.hh index 83bff63f704cdb0eef6a58a4dfdceb2f637cddf5..793dea331e794f10e02021d2e6d7c90832c92c6d 100644 --- a/dumux/discretization/cellcentered/tpfa/computetransmissibility.hh +++ b/dumux/discretization/cellcentered/tpfa/computetransmissibility.hh @@ -44,6 +44,7 @@ namespace Dumux { * \param extrusionFactor The extrusion factor of the scv */ template< class SubControlVolumeFace, class SubControlVolume, class Tensor > +[[deprecated("Will be removed after release 3.6. Use interface with additional fvGeometry parameter instead.")]] typename Tensor::field_type computeTpfaTransmissibility(const SubControlVolumeFace& scvf, const SubControlVolume& scv, const Tensor& T, @@ -60,6 +61,37 @@ typename Tensor::field_type computeTpfaTransmissibility(const SubControlVolumeFa return (Knormal*distanceVector) * extrusionFactor; } +/*! + * \ingroup CCTpfaDiscretization + * \brief Free function to evaluate the Tpfa transmissibility + * associated with the flux (in the form of flux = T*gradU) across a + * sub-control volume face stemming from a given sub-control + * volume with corresponding tensor T. + * + * \param fvGeometry The element-centered control volume geometry + * \param scvf The sub-control volume face + * \param scv The neighboring sub-control volume + * \param T The tensor living in the neighboring scv + * \param extrusionFactor The extrusion factor of the scv + */ +template< class FVElementGeometry, class Tensor > +typename Tensor::field_type computeTpfaTransmissibility(const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const typename FVElementGeometry::SubControlVolume& scv, + const Tensor& T, + typename FVElementGeometry::SubControlVolume::Traits::Scalar extrusionFactor) +{ + using GlobalPosition = typename FVElementGeometry::SubControlVolumeFace::Traits::GlobalPosition; + GlobalPosition Knormal; + T.mv(scvf.unitOuterNormal(), Knormal); + + auto distanceVector = scvf.ipGlobal(); + distanceVector -= scv.center(); + distanceVector /= distanceVector.two_norm2(); + + return (Knormal*distanceVector) * extrusionFactor; +} + /*! * \ingroup CCTpfaDiscretization * \brief Free function to evaluate the Tpfa transmissibility @@ -76,6 +108,7 @@ template< class SubControlVolumeFace, class SubControlVolume, class Tensor, typename std::enable_if_t::value, int> = 0 > +[[deprecated("Will be removed after release 3.6. Use interface with additional fvGeometry parameter instead.")]] Tensor computeTpfaTransmissibility(const SubControlVolumeFace& scvf, const SubControlVolume &scv, Tensor t, @@ -88,6 +121,35 @@ Tensor computeTpfaTransmissibility(const SubControlVolumeFace& scvf, return t * extrusionFactor * (distanceVector * scvf.unitOuterNormal()); } +/*! + * \ingroup CCTpfaDiscretization + * \brief Free function to evaluate the Tpfa transmissibility + * associated with the flux (in the form of flux = T*gradU) across a + * sub-control volume face stemming from a given sub-control + * volume for the case where T is just a scalar + * + * \param fvGeometry The element-centered control volume geometry + * \param scvf The sub-control volume face + * \param scv The neighboring sub-control volume + * \param t The scalar quantity living in the neighboring scv + * \param extrusionFactor The extrusion factor of the scv + */ +template< class FVElementGeometry, + class Tensor, + typename std::enable_if_t::value, int> = 0 > +Tensor computeTpfaTransmissibility(const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const typename FVElementGeometry::SubControlVolume& scv, + Tensor t, + typename FVElementGeometry::SubControlVolumeFace::Traits::Scalar extrusionFactor) +{ + auto distanceVector = scvf.ipGlobal(); + distanceVector -= scv.center(); + distanceVector /= distanceVector.two_norm2(); + + return t * extrusionFactor * (distanceVector * scvf.unitOuterNormal()); +} + } // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh index 6129f321aa8db67f55b27babd75fbfb0fe7d9a91..e36175f2f4e3f82b969f456d4a117d4a7746ca5a 100644 --- a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh +++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh @@ -194,6 +194,18 @@ public: bool hasBoundaryScvf() const { return gridGeometry().hasBoundaryScvf(scvIndices_[0]); } + typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const + { + assert(isBound()); + return scv.geometry(); + } + + typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const + { + assert(isBound()); + return scvf.geometry(); + } + private: std::optional element_; @@ -357,6 +369,18 @@ public: bool hasBoundaryScvf() const { return hasBoundaryScvf_; } + typename SubControlVolume::Traits::Geometry geometry (const SubControlVolume& scv) const + { + assert(isBound()); + return scv.geometry(); + } + + typename SubControlVolumeFace::Traits::Geometry geometry (const SubControlVolumeFace& scvf) const + { + assert(isBound()); + return scvf.geometry(); + } + private: //! Binding of an element preparing the geometries of the whole stencil //! called by the local jacobian to prepare element assembly diff --git a/dumux/discretization/extrusion.hh b/dumux/discretization/extrusion.hh index 19fabb9e85f63cef6655480a9b11223ce21c52e7..aa30f7ffbd8c4d964495cfdcf267cc97e2132923 100644 --- a/dumux/discretization/extrusion.hh +++ b/dumux/discretization/extrusion.hh @@ -48,12 +48,12 @@ struct NoExtrusion static constexpr auto volume(const SCV& scv) { return scv.volume(); } - template - static constexpr auto area(const FVGeo& fvGeometry, const typename FVGeo::SubControlVolumeFace& scvf) + template + static constexpr auto area(const FVGeo&, const SCVF& scvf) { return scvf.area(); } - template - static constexpr auto volume(const FVGeo& fvGeometry, const typename FVGeo::SubControlVolume& scv) + template + static constexpr auto volume(const FVGeo&, const SCV& scv) { return scv.volume(); } template @@ -87,10 +87,9 @@ struct RotationalExtrusion return scvf.area()*2.0*M_PI*scvf.center()[radialAxis]; } - template - static constexpr auto area(const FVGeo& fvGeometry, const typename FVGeo::SubControlVolumeFace& scvf) + template + static constexpr auto area(const FVGeo&, const SCVF& scvf) { - using SCVF = typename FVGeo::SubControlVolumeFace; static_assert(int(SCVF::Traits::Geometry::mydimension) == int(SCVF::Traits::Geometry::coorddimension-1), "Area element to be called with a codim-1-entity!"); static_assert(SCVF::Traits::Geometry::coorddimension <= 2, "Axis rotation only makes sense for geometries up to 2D!"); static_assert(radialAxis < int(SCVF::Traits::Geometry::coorddimension), "Illegal radial axis!"); @@ -115,10 +114,9 @@ struct RotationalExtrusion return scv.volume()*2.0*M_PI*scv.center()[radialAxis]; } - template - static constexpr auto volume(const FVGeo& fvGeometry, const typename FVGeo::SubControlVolume& scv) + template + static constexpr auto volume(const FVGeo&, const SCV& scv) { - using SCV = typename FVGeo::SubControlVolume; static_assert(int(SCV::Traits::Geometry::mydimension) == int(SCV::Traits::Geometry::coorddimension), "Volume element to be called with a codim-0-entity!"); static_assert(SCV::Traits::Geometry::coorddimension <= 2, "Axis rotation only makes sense for geometries up to 2D!"); static_assert(radialAxis < int(SCV::Traits::Geometry::coorddimension), "Illegal radial axis!"); @@ -163,10 +161,9 @@ struct SphericalExtrusion return 4.0*M_PI*radius*radius; } - template - static constexpr auto area(const FVGeo& fvGeometry, const typename FVGeo::SubControlVolumeFace& scvf) + template + static constexpr auto area(const FVGeo&, const SCVF& scvf) { - using SCVF = typename FVGeo::SubControlVolumeFace; static_assert(int(SCVF::Traits::Geometry::mydimension) == int(SCVF::Traits::Geometry::coorddimension-1), "Area element to be called with a codim-1-entity!"); static_assert(SCVF::Traits::Geometry::coorddimension == 1, "Spherical rotation only makes sense for 1D geometries!"); @@ -193,10 +190,9 @@ struct SphericalExtrusion return 4.0/3.0*M_PI*abs(radius1*radius1*radius1 - radius0*radius0*radius0); } - template - static constexpr auto volume(const FVGeo& fvGeometry, const typename FVGeo::SubControlVolume& scv) + template + static constexpr auto volume(const FVGeo& fvGeometry, const SCV& scv) { - using SCV = typename FVGeo::SubControlVolume; static_assert(int(SCV::Traits::Geometry::mydimension) == int(SCV::Traits::Geometry::coorddimension), "Volume element to be called with a codim-0-entity!"); static_assert(SCV::Traits::Geometry::coorddimension == 1, "Spherical rotation only makes sense for 1D geometries!"); diff --git a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh index 9f08c7e6b17cf802fd719f82add21f90e8f33052..2f62de3574fb52b6f52b3ef23f19c4368a014c96 100644 --- a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh +++ b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh @@ -221,6 +221,12 @@ public: DUNE_THROW(Dune::InvalidStateException, "No outside scvf found"); } + typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const + { return scv.geometry(); } + + typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const + { return scvf.geometry(); } + private: const auto& scvfIndices_() const diff --git a/dumux/discretization/staggered/fvelementgeometry.hh b/dumux/discretization/staggered/fvelementgeometry.hh index 6581a68b9c8e264289fc3638357247702e6bf493..0f37a115503b188e54a7b27ba945fffca03d0dac 100644 --- a/dumux/discretization/staggered/fvelementgeometry.hh +++ b/dumux/discretization/staggered/fvelementgeometry.hh @@ -84,7 +84,7 @@ public: * \brief Base class for the finite volume geometry vector for staggered models * This locally builds up the sub control volumes and sub control volume faces * for each element. - * Specialization for grid caching enabled + * Specialization for grid caching disabled * \tparam GG the finite volume grid geometry type */ template @@ -219,6 +219,12 @@ public: bool hasBoundaryScvf() const { return hasBoundaryScvf_; } + typename SubControlVolumeFace::Traits::Geometry geometry (const SubControlVolumeFace& scvf) const + { + assert(isBound()); + return scvf.geometry(); + } + private: //! Binding of an element preparing the geometries only inside the element void bindElement_(const Element& element) diff --git a/dumux/discretization/walldistance.hh b/dumux/discretization/walldistance.hh index 045b8e2ce6e6df39fb1be1da054007227c361367..ae485bee92b3d7ef7df1610b617b2c1d7c2c6f63 100644 --- a/dumux/discretization/walldistance.hh +++ b/dumux/discretization/walldistance.hh @@ -197,7 +197,7 @@ private: { if (scvf.boundary() && considerFace(fvGeometry, scvf)) { - const auto& geo = scvf.geometry(); + const auto& geo = fvGeometry.geometry(scvf); CornerStorage corners; for (int i = 0; i < geo.corners(); ++i) corners.push_back(geo.corner(i)); diff --git a/dumux/flux/box/darcyslaw.hh b/dumux/flux/box/darcyslaw.hh index e01000f1b39f0ce27c0a0b5a5e92df2e617114ab..ceb511b949740f6be254ce4a82d4e256fd26ec52 100644 --- a/dumux/flux/box/darcyslaw.hh +++ b/dumux/flux/box/darcyslaw.hh @@ -131,7 +131,7 @@ public: gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center())); // apply the permeability and return the flux - return -1.0*vtmv(scvf.unitOuterNormal(), K, gradP)*Extrusion::area(scvf); + return -1.0*vtmv(scvf.unitOuterNormal(), K, gradP)*Extrusion::area(fvGeometry, scvf); } // compute transmissibilities ti for analytical jacobians @@ -160,7 +160,7 @@ public: std::vector ti(fvGeometry.numScv()); for (const auto& scv : scvs(fvGeometry)) ti[scv.indexInElement()] = - -1.0*Extrusion::area(scvf)*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.indexInElement())); + -1.0*Extrusion::area(fvGeometry, scvf)*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.indexInElement())); return ti; } diff --git a/dumux/flux/box/dispersionflux.hh b/dumux/flux/box/dispersionflux.hh index b1a6c6d1cac254802d17d3e2133c344d20f8196b..0a3b7dfe81c1888725250c5a83eb3a4454cb5926 100644 --- a/dumux/flux/box/dispersionflux.hh +++ b/dumux/flux/box/dispersionflux.hh @@ -164,7 +164,7 @@ public: gradTemp.axpy(elemVolVars[scv].temperature(), fluxVarsCache.gradN(scv.indexInElement())); // compute the heat conduction flux - return -1.0*vtmv(scvf.unitOuterNormal(), dispersionTensor, gradTemp)*Extrusion::area(scvf); + return -1.0*vtmv(scvf.unitOuterNormal(), dispersionTensor, gradTemp)*Extrusion::area(fvGeometry, scvf); } }; diff --git a/dumux/flux/box/effectivestresslaw.hh b/dumux/flux/box/effectivestresslaw.hh index ebd8d4ac594d5530fc1a9bca2c0a9709e230209c..de8271e9e4e0d7f685b1b4844da518d291459b0c 100644 --- a/dumux/flux/box/effectivestresslaw.hh +++ b/dumux/flux/box/effectivestresslaw.hh @@ -81,7 +81,7 @@ public: ForceVector scvfForce(0.0); sigma.mv(scvf.unitOuterNormal(), scvfForce); - scvfForce *= Extrusion::area(scvf); + scvfForce *= Extrusion::area(fvGeometry, scvf); return scvfForce; } diff --git a/dumux/flux/box/fickslaw.hh b/dumux/flux/box/fickslaw.hh index 46c5badcf83d0049272fc98d3273c48d38105b59..f94549f3ee2bc6049923d09c3c8406186635a31f 100644 --- a/dumux/flux/box/fickslaw.hh +++ b/dumux/flux/box/fickslaw.hh @@ -162,7 +162,7 @@ public: ti[compIdx].resize(fvGeometry.numScv()); for (auto&& scv : scvs(fvGeometry)) - ti[compIdx][scv.indexInElement()] = -rho*vtmv(scvf.unitOuterNormal(), diffCoeff, fluxVarCache.gradN(scv.indexInElement()))*Extrusion::area(scvf); + ti[compIdx][scv.indexInElement()] = -rho*vtmv(scvf.unitOuterNormal(), diffCoeff, fluxVarCache.gradN(scv.indexInElement()))*Extrusion::area(fvGeometry, scvf); } return ti; @@ -214,7 +214,7 @@ private: Dune::FieldVector gradX(0.0); for (auto&& scv : scvs(fvGeometry)) gradX.axpy(massOrMoleFraction(scv), fluxVarsCache.gradN(scv.indexInElement())); - return -1.0*preFactor*vtmv(scvf.unitOuterNormal(), D, gradX)*Extrusion::area(scvf); + return -1.0*preFactor*vtmv(scvf.unitOuterNormal(), D, gradX)*Extrusion::area(fvGeometry, scvf); } }; diff --git a/dumux/flux/box/forchheimerslaw.hh b/dumux/flux/box/forchheimerslaw.hh index f7c128864257c3d12af6fa40d410bc6ad3e9a89d..ae6e3decb11886b37abb8cfd878eedc5a4daf253 100644 --- a/dumux/flux/box/forchheimerslaw.hh +++ b/dumux/flux/box/forchheimerslaw.hh @@ -158,7 +158,7 @@ public: darcyVelocity); Scalar flux = velocity * scvf.unitOuterNormal(); - flux *= Extrusion::area(scvf); + flux *= Extrusion::area(fvGeometry, scvf); return flux; } diff --git a/dumux/flux/box/fourierslaw.hh b/dumux/flux/box/fourierslaw.hh index e91d6503a35af712c51000b88ef1fee140dba760..f29d3b946060407108fec9b2055f3d7dca506d69 100644 --- a/dumux/flux/box/fourierslaw.hh +++ b/dumux/flux/box/fourierslaw.hh @@ -96,7 +96,7 @@ public: gradTemp.axpy(elemVolVars[scv].temperature(), fluxVarsCache.gradN(scv.indexInElement())); // compute the heat conduction flux - return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*Extrusion::area(scvf); + return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*Extrusion::area(fvGeometry, scvf); } }; diff --git a/dumux/flux/box/fourierslawnonequilibrium.hh b/dumux/flux/box/fourierslawnonequilibrium.hh index 1c5024088b1a29ed91c731db106e1bd83e7b5f36..f0118f6b8298cd6e070fb6ddafdb569a61c2ecd0 100644 --- a/dumux/flux/box/fourierslawnonequilibrium.hh +++ b/dumux/flux/box/fourierslawnonequilibrium.hh @@ -119,7 +119,7 @@ public: } // compute the heat conduction flux - return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*Extrusion::area(scvf); + return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*Extrusion::area(fvGeometry, scvf); } }; diff --git a/dumux/flux/box/hookeslaw.hh b/dumux/flux/box/hookeslaw.hh index f750acd26f406ce194ac6f516d49a8773decaeea..d29dfe167d88cdcf7bdb8cef642aeac4e592491d 100644 --- a/dumux/flux/box/hookeslaw.hh +++ b/dumux/flux/box/hookeslaw.hh @@ -81,7 +81,7 @@ public: ForceVector scvfForce(0.0); sigma.mv(scvf.unitOuterNormal(), scvfForce); - scvfForce *= Extrusion::area(scvf); + scvfForce *= Extrusion::area(fvGeometry, scvf); return scvfForce; } diff --git a/dumux/flux/box/maxwellstefanslaw.hh b/dumux/flux/box/maxwellstefanslaw.hh index 4e21b7eb1750dc8f28e23b0dc84f355c33801ac7..cdea991878ed89b9a9f4e9f6b810dd9a00da85ff 100644 --- a/dumux/flux/box/maxwellstefanslaw.hh +++ b/dumux/flux/box/maxwellstefanslaw.hh @@ -138,7 +138,7 @@ public: normalX[compIdx] = gradX *scvf.unitOuterNormal(); } reducedDiffusionMatrix.solve(reducedFlux,normalX); - reducedFlux *= -1.0*rho*Extrusion::area(scvf); + reducedFlux *= -1.0*rho*Extrusion::area(fvGeometry, scvf); for (int compIdx = 0; compIdx < numComponents-1; compIdx++) { diff --git a/dumux/flux/cctpfa/darcyslaw.hh b/dumux/flux/cctpfa/darcyslaw.hh index fa9b7eb164d3e07f112da1952d70c41352ff40bd..378e79b6745bcdebde5dd963c1bb104eda0e4e69 100644 --- a/dumux/flux/cctpfa/darcyslaw.hh +++ b/dumux/flux/cctpfa/darcyslaw.hh @@ -196,7 +196,7 @@ class CCTpfaDarcysLaw //! compute alpha := n^T*K*g const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor(); - Scalar flux = tij*(pInside - pOutside) + rho*Extrusion::area(scvf)*alpha_inside; + Scalar flux = tij*(pInside - pOutside) + rho*Extrusion::area(fvGeometry, scvf)*alpha_inside; //! On interior faces we have to add K-weighted gravitational contributions if (!scvf.boundary()) @@ -204,8 +204,8 @@ class CCTpfaDarcysLaw const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); const auto outsideK = outsideVolVars.permeability(); const auto outsideTi = fvGeometry.gridGeometry().isPeriodic() - ? computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, outsideK, outsideVolVars.extrusionFactor()) - : -1.0*computeTpfaTransmissibility(scvf, outsideScv, outsideK, outsideVolVars.extrusionFactor()); + ? computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()), outsideScv, outsideK, outsideVolVars.extrusionFactor()) + : -1.0*computeTpfaTransmissibility(fvGeometry, scvf, outsideScv, outsideK, outsideVolVars.extrusionFactor()); const auto alpha_outside = vtmv(scvf.unitOuterNormal(), outsideK, g)*outsideVolVars.extrusionFactor(); flux -= rho*tij/outsideTi*(alpha_inside - alpha_outside); @@ -239,13 +239,13 @@ class CCTpfaDarcysLaw const auto& insideScv = fvGeometry.scv(insideScvIdx); const auto& insideVolVars = elemVolVars[insideScvIdx]; - const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, + const Scalar ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, getPermeability_(problem, insideVolVars, scvf.ipGlobal()), insideVolVars.extrusionFactor()); // on the boundary (dirichlet) we only need ti if (scvf.boundary()) - tij = Extrusion::area(scvf)*ti; + tij = Extrusion::area(fvGeometry, scvf)*ti; // otherwise we compute a tpfa harmonic mean else @@ -256,15 +256,15 @@ class CCTpfaDarcysLaw const auto& outsideScv = fvGeometry.scv(outsideScvIdx); const auto& outsideVolVars = elemVolVars[outsideScvIdx]; const Scalar tj = fvGeometry.gridGeometry().isPeriodic() - ? computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor()) - : -1.0*computeTpfaTransmissibility(scvf, outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor()); + ? computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()), outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor()) + : -1.0*computeTpfaTransmissibility(fvGeometry, scvf, outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor()); // harmonic mean (check for division by zero!) // TODO: This could lead to problems!? Is there a better way to do this? if (ti*tj <= 0.0) tij = 0; else - tij = Extrusion::area(scvf)*(ti * tj)/(ti + tj); + tij = Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj); } return tij; @@ -390,7 +390,7 @@ public: Scalar sumPTi(tij*pInside); // add inside gravitational contribution - sumPTi += rho*Extrusion::area(scvf) + sumPTi += rho*Extrusion::area(fvGeometry, scvf) *insideVolVars.extrusionFactor() *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g); @@ -404,7 +404,7 @@ public: sumPTi += outsideFluxVarsCache.advectionTij()*outsideVolVars.pressure(phaseIdx); // add outside gravitational contribution - sumPTi += rho*Extrusion::area(scvf) + sumPTi += rho*Extrusion::area(fvGeometry, scvf) *outsideVolVars.extrusionFactor() *vtmv(flippedScvf.unitOuterNormal(), outsideVolVars.permeability(), g); } @@ -415,7 +415,7 @@ public: //! precompute alpha := n^T*K*g const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor(); - Scalar flux = tij*(pInside - pOutside) + Extrusion::area(scvf)*rho*alpha_inside; + Scalar flux = tij*(pInside - pOutside) + Extrusion::area(fvGeometry, scvf)*rho*alpha_inside; //! On interior faces with one neighbor we have to add K-weighted gravitational contributions if (!scvf.boundary() && scvf.numOutsideScvs() == 1) @@ -423,7 +423,7 @@ public: const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); const auto& outsideScvf = fvGeometry.flipScvf(scvf.index()); const auto outsideK = outsideVolVars.permeability(); - const auto outsideTi = computeTpfaTransmissibility(outsideScvf, outsideScv, outsideK, outsideVolVars.extrusionFactor()); + const auto outsideTi = computeTpfaTransmissibility(fvGeometry, outsideScvf, outsideScv, outsideK, outsideVolVars.extrusionFactor()); const auto alpha_outside = vtmv(outsideScvf.unitOuterNormal(), outsideK, g)*outsideVolVars.extrusionFactor(); flux -= rho*tij/outsideTi*(alpha_inside + alpha_outside); @@ -481,13 +481,13 @@ public: const auto& insideScv = fvGeometry.scv(insideScvIdx); const auto& insideVolVars = elemVolVars[insideScvIdx]; - const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, + const Scalar ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, getPermeability_(problem, insideVolVars, scvf.ipGlobal()), insideVolVars.extrusionFactor()); // for the boundary (dirichlet) or at branching points we only need ti if (scvf.boundary() || scvf.numOutsideScvs() > 1) - tij = Extrusion::area(scvf)*ti; + tij = Extrusion::area(fvGeometry, scvf)*ti; // otherwise we compute a tpfa harmonic mean else @@ -497,7 +497,7 @@ public: // refers to the scv of our element, so we use the scv method const auto& outsideScv = fvGeometry.scv(outsideScvIdx); const auto& outsideVolVars = elemVolVars[outsideScvIdx]; - const Scalar tj = computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, + const Scalar tj = computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()), outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor()); @@ -506,7 +506,7 @@ public: if (ti*tj <= 0.0) tij = 0; else - tij = Extrusion::area(scvf)*(ti * tj)/(ti + tj); + tij = Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj); } return tij; diff --git a/dumux/flux/cctpfa/dispersionflux.hh b/dumux/flux/cctpfa/dispersionflux.hh index e3428afa204e2113c185abf7d25cfcc2aaeef64c..9c0fe1e5d2506ee8c5501fc7bf09197d0e915354 100644 --- a/dumux/flux/cctpfa/dispersionflux.hh +++ b/dumux/flux/cctpfa/dispersionflux.hh @@ -123,7 +123,7 @@ public: ModelTraits::CompositionalDispersionModel::compositionalDispersionTensor(problem, scvf, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx, compIdx); - const auto dij = computeTpfaTransmissibility(scvf, fvGeometry.scv(scvf.insideScvIdx()), dispersionTensor, insideVolVars.extrusionFactor()); + const auto dij = computeTpfaTransmissibility(fvGeometry, scvf, fvGeometry.scv(scvf.insideScvIdx()), dispersionTensor, insideVolVars.extrusionFactor()); const auto xInside = massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx); const auto xOutide = massOrMoleFraction(outsideVolVars, referenceSystem, phaseIdx, compIdx); diff --git a/dumux/flux/cctpfa/fickslaw.hh b/dumux/flux/cctpfa/fickslaw.hh index ce5847b4b83dc8a75e6b6e3ed1951f3e4d52a2c4..7d681adddbf66d982b79d8a18a1b7a22262839e5 100644 --- a/dumux/flux/cctpfa/fickslaw.hh +++ b/dumux/flux/cctpfa/fickslaw.hh @@ -241,12 +241,12 @@ public: const auto insideDiffCoeff = getDiffCoeff(insideVolVars); - const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, insideDiffCoeff, insideVolVars.extrusionFactor()); + const Scalar ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, insideDiffCoeff, insideVolVars.extrusionFactor()); // for the boundary (dirichlet) or at branching points we only need ti Scalar tij; if (scvf.boundary() || scvf.numOutsideScvs() > 1) - tij = Extrusion::area(scvf)*ti; + tij = Extrusion::area(fvGeometry, scvf)*ti; // otherwise we compute a tpfa harmonic mean else @@ -259,9 +259,9 @@ public: Scalar tj; if constexpr (dim == dimWorld) // assume the normal vector from outside is anti parallel so we save flipping a vector - tj = -1.0*computeTpfaTransmissibility(scvf, outsideScv, outsideDiffCoeff, outsideVolVars.extrusionFactor()); + tj = -1.0*computeTpfaTransmissibility(fvGeometry, scvf, outsideScv, outsideDiffCoeff, outsideVolVars.extrusionFactor()); else - tj = computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), + tj = computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()), outsideScv, outsideDiffCoeff, outsideVolVars.extrusionFactor()); @@ -270,7 +270,7 @@ public: if (ti*tj <= 0.0) tij = 0; else - tij = Extrusion::area(scvf)*(ti * tj)/(ti + tj); + tij = Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj); } return tij; diff --git a/dumux/flux/cctpfa/forchheimerslaw.hh b/dumux/flux/cctpfa/forchheimerslaw.hh index 7e1065254d36b73fb6d62bda8c381e74bb6a1a8b..5e05afe757dcd185b6fce760346ab777261f9453 100644 --- a/dumux/flux/cctpfa/forchheimerslaw.hh +++ b/dumux/flux/cctpfa/forchheimerslaw.hh @@ -186,7 +186,7 @@ class CCTpfaForchheimersLaw 1) { - tij = Extrusion::area(scvf)*ti; + tij = Extrusion::area(fvGeometry, scvf)*ti; } // otherwise we compute a tpfa harmonic mean else @@ -201,15 +201,15 @@ public: Scalar tj; if constexpr (dim == dimWorld) // assume the normal vector from outside is anti parallel so we save flipping a vector - tj = -1.0*computeTpfaTransmissibility(scvf, outsideScv, outsideLambda, outsideVolVars.extrusionFactor()); + tj = -1.0*computeTpfaTransmissibility(fvGeometry, scvf, outsideScv, outsideLambda, outsideVolVars.extrusionFactor()); else - tj = computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, outsideLambda, outsideVolVars.extrusionFactor()); + tj = computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()), outsideScv, outsideLambda, outsideVolVars.extrusionFactor()); // check for division by zero! if (ti*tj <= 0.0) tij = 0; else - tij = Extrusion::area(scvf)*(ti * tj)/(ti + tj); + tij = Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj); } return tij; diff --git a/dumux/flux/cctpfa/fourierslawnonequilibrium.hh b/dumux/flux/cctpfa/fourierslawnonequilibrium.hh index a2ddb1e83d8db26a494b14825936567240d26660..d92e1333d290ab0a32ca662bc8b22dd72878f835 100644 --- a/dumux/flux/cctpfa/fourierslawnonequilibrium.hh +++ b/dumux/flux/cctpfa/fourierslawnonequilibrium.hh @@ -125,11 +125,11 @@ public: }; const auto insideLambda = computeLambda(insideVolVars); - const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, insideLambda, insideVolVars.extrusionFactor()); + const Scalar ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, insideLambda, insideVolVars.extrusionFactor()); // for the boundary (dirichlet) or at branching points we only need ti if (scvf.boundary() || scvf.numOutsideScvs() > 1) - return Extrusion::area(scvf)*ti; + return Extrusion::area(fvGeometry, scvf)*ti; else // otherwise we compute a tpfa harmonic mean { const auto outsideScvIdx = scvf.outsideScvIdx(); @@ -140,15 +140,15 @@ public: Scalar tj; if (dim == dimWorld) // assume the normal vector from outside is anti parallel so we save flipping a vector - tj = -1.0*computeTpfaTransmissibility(scvf, outsideScv, outsideLambda, outsideVolVars.extrusionFactor()); + tj = -1.0*computeTpfaTransmissibility(fvGeometry, scvf, outsideScv, outsideLambda, outsideVolVars.extrusionFactor()); else - tj = computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, outsideLambda, outsideVolVars.extrusionFactor()); + tj = computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()), outsideScv, outsideLambda, outsideVolVars.extrusionFactor()); // check for division by zero! if (ti*tj <= 0.0) return 0.0; else - return Extrusion::area(scvf)*(ti * tj)/(ti + tj); + return Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj); } } }; diff --git a/dumux/flux/cctpfa/maxwellstefanslaw.hh b/dumux/flux/cctpfa/maxwellstefanslaw.hh index c45578cb0c7f0468253c9f5a75aae5c7fc10018b..9e08cb09eb886baaa142cc90555a16ac89e45ac4 100644 --- a/dumux/flux/cctpfa/maxwellstefanslaw.hh +++ b/dumux/flux/cctpfa/maxwellstefanslaw.hh @@ -195,7 +195,7 @@ public: reducedDiffusionMatrixInside.mv(helperVector, reducedFlux); } - reducedFlux *= -Extrusion::area(scvf); + reducedFlux *= -Extrusion::area(fvGeometry, scvf); for (int compIdx = 0; compIdx < numComponents-1; compIdx++) { componentFlux[compIdx] = reducedFlux[compIdx]; diff --git a/dumux/flux/cvfe/darcyslaw.hh b/dumux/flux/cvfe/darcyslaw.hh index 42ec2a1047e716707dba855e6318cbe19a8eca32..ed05f5c2973d7e50b13883b5f866081ddd3e7304 100644 --- a/dumux/flux/cvfe/darcyslaw.hh +++ b/dumux/flux/cvfe/darcyslaw.hh @@ -112,7 +112,7 @@ public: gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center())); // apply the permeability and return the flux - return -1.0*vtmv(scvf.unitOuterNormal(), K, gradP)*Extrusion::area(scvf); + return -1.0*vtmv(scvf.unitOuterNormal(), K, gradP)*Extrusion::area(fvGeometry, scvf); } // compute transmissibilities ti for analytical Jacobians @@ -141,7 +141,7 @@ public: std::vector ti(fvGeometry.numScv()); for (const auto& scv : scvs(fvGeometry)) ti[scv.indexInElement()] = - -1.0*Extrusion::area(scvf)*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.indexInElement())); + -1.0*Extrusion::area(fvGeometry, scvf)*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.indexInElement())); return ti; } diff --git a/dumux/flux/staggered/freeflow/fickslaw.hh b/dumux/flux/staggered/freeflow/fickslaw.hh index 96359c313dff2711e6e3f992b04bd2a96b7a0db4..ae4a8b1885590a25ef1f5a046cdb2beef0fe2c17 100644 --- a/dumux/flux/staggered/freeflow/fickslaw.hh +++ b/dumux/flux/staggered/freeflow/fickslaw.hh @@ -150,7 +150,7 @@ public: const Scalar cumulativeFlux = std::accumulate(flux.begin(), flux.end(), 0.0); flux[FluidSystem::getMainComponent(0)] = -cumulativeFlux; - flux *= Extrusion::area(scvf); + flux *= Extrusion::area(fvGeometry, scvf); return flux; } diff --git a/dumux/flux/staggered/freeflow/fourierslaw.hh b/dumux/flux/staggered/freeflow/fourierslaw.hh index 1221ec07c1674065d305e775aa9e53da13ed41b6..76d2e2230c9f27175472b9a1610f671041b341a3 100644 --- a/dumux/flux/staggered/freeflow/fourierslaw.hh +++ b/dumux/flux/staggered/freeflow/fourierslaw.hh @@ -106,7 +106,7 @@ public: flux = avgLambda * (insideTemperature - outsideTemperature) / (insideDistance + outsideDistance); } - flux *= Extrusion::area(scvf); + flux *= Extrusion::area(fvGeometry, scvf); return flux; } }; diff --git a/dumux/flux/staggered/freeflow/maxwellstefanslaw.hh b/dumux/flux/staggered/freeflow/maxwellstefanslaw.hh index d6408f78930a02e9d68cfe80a23ef9175446d173..66960d446eb26d8d12e2a964508fdb14beee3648 100644 --- a/dumux/flux/staggered/freeflow/maxwellstefanslaw.hh +++ b/dumux/flux/staggered/freeflow/maxwellstefanslaw.hh @@ -196,7 +196,7 @@ public: reducedDiffusionMatrixInside.mv(helperVector, reducedFlux); } - reducedFlux *= -Extrusion::area(scvf); + reducedFlux *= -Extrusion::area(fvGeometry, scvf); for (int compIdx = 0; compIdx < numComponents-1; compIdx++) { diff --git a/dumux/freeflow/compositional/staggered/fluxvariables.hh b/dumux/freeflow/compositional/staggered/fluxvariables.hh index b6c3db8f517e7cb53177af4067e5d4d46e4c47dd..cf9bf0659ff3bf0381a0d97232be7de74db5377a 100644 --- a/dumux/freeflow/compositional/staggered/fluxvariables.hh +++ b/dumux/freeflow/compositional/staggered/fluxvariables.hh @@ -88,7 +88,7 @@ public: return density * fraction; }; - flux[compIdx] = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTerm); + flux[compIdx] = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTerm); //check for the reference system and adapt units of the diffusive flux accordingly. if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged) diff --git a/dumux/freeflow/navierstokes/momentum/diamond/flux.hh b/dumux/freeflow/navierstokes/momentum/diamond/flux.hh index e1b47b96f7ca2ddabee4fdc6d5398c018c485784..03c3421c9f8fb0e0cc7749d488fb1b46ddf6fe88 100644 --- a/dumux/freeflow/navierstokes/momentum/diamond/flux.hh +++ b/dumux/freeflow/navierstokes/momentum/diamond/flux.hh @@ -145,7 +145,7 @@ public: static const auto upwindWeight = getParamFromGroup(context.problem().paramGroup(), "Flux.UpwindWeight"); const auto advectiveTermIntegrand = density*vn * (upwindWeight * upwindVelocity + (1.0-upwindWeight)*downwindVelocity); - return advectiveTermIntegrand * Extrusion::area(scvf) * insideVolVars.extrusionFactor(); + return advectiveTermIntegrand * Extrusion::area(fvGeometry, scvf) * insideVolVars.extrusionFactor(); } /*! @@ -186,7 +186,7 @@ public: if (enableDilatationTerm) diffusiveFlux += 2.0/3.0 * mu * trace(gradV) * scvf.unitOuterNormal(); - diffusiveFlux *= Extrusion::area(scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor(); + diffusiveFlux *= Extrusion::area(fvGeometry, scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor(); return diffusiveFlux; } @@ -208,7 +208,7 @@ public: const auto referencePressure = context.problem().referencePressure(); NumEqVector pn(scvf.unitOuterNormal()); - pn *= (pressure-referencePressure)*Extrusion::area(scvf)*elemVolVars[scvf.insideScvIdx()].extrusionFactor(); + pn *= (pressure-referencePressure)*Extrusion::area(fvGeometry, scvf)*elemVolVars[scvf.insideScvIdx()].extrusionFactor(); return pn; } diff --git a/dumux/freeflow/navierstokes/momentum/fluxvariables.hh b/dumux/freeflow/navierstokes/momentum/fluxvariables.hh index 4d350f9e3e3355107f2528226cf0fa6cac68c687..6cc1f97c56a88f46eda33293dd28cf02fa3bc06a 100644 --- a/dumux/freeflow/navierstokes/momentum/fluxvariables.hh +++ b/dumux/freeflow/navierstokes/momentum/fluxvariables.hh @@ -185,7 +185,7 @@ public: static const Scalar factor = enableUnsymmetrizedVelocityGradient ? 1.0 : 2.0; const auto mu = this->problem().effectiveViscosity(this->element(), this->fvGeometry(), this->scvFace()); - result -= factor * mu * velGradII * Extrusion::area(scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor() * scvf.directionSign(); + result -= factor * mu * velGradII * Extrusion::area(fvGeometry, scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor() * scvf.directionSign(); static const bool enableDilatationTerm = getParamFromGroup(this->problem().paramGroup(), "FreeFlow.EnableDilatationTerm", false); if (enableDilatationTerm) @@ -199,7 +199,7 @@ public: divergence += VelocityGradients::velocityGradII(fvGeometry, otherFrontalScvf, elemVolVars); } - result += 2.0/3.0 * mu * divergence * scvf.directionSign() * Extrusion::area(scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor(); + result += 2.0/3.0 * mu * divergence * scvf.directionSign() * Extrusion::area(fvGeometry, scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor(); } return result; @@ -262,7 +262,7 @@ public: } // Account for the area of the staggered lateral face. - return result * Extrusion::area(scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor(); + return result * Extrusion::area(fvGeometry, scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor(); } /*! @@ -290,7 +290,7 @@ public: // The pressure force needs to take the extruded scvf area into account. const auto pressure = this->problem().pressure(this->element(), this->fvGeometry(), scvf); - result = pressure*Extrusion::area(scvf)*this->elemVolVars()[scvf.insideScvIdx()].extrusionFactor(); + result = pressure*Extrusion::area(this->fvGeometry(), scvf)*this->elemVolVars()[scvf.insideScvIdx()].extrusionFactor(); // The pressure contribution calculated above might have a much larger numerical value compared to the viscous or inertial forces. // This may lead to numerical inaccuracies due to loss of significance (cancellantion) for the final residual value. @@ -347,7 +347,7 @@ public: const Scalar transportedMomentum = selfIsUpstream ? (upwindWeight * velocitySelf + (1.0 - upwindWeight) * velocityOpposite) * density : (upwindWeight * velocityOpposite + (1.0 - upwindWeight) * velocitySelf) * density; - return transportingVelocity * transportedMomentum * scvf.directionSign() * Extrusion::area(scvf) * extrusionFactor_(elemVolVars, scvf); + return transportingVelocity * transportedMomentum * scvf.directionSign() * Extrusion::area(this->fvGeometry(), scvf) * extrusionFactor_(elemVolVars, scvf); } /*! @@ -470,7 +470,7 @@ public: : (upwindWeight * outsideMomentum + (1.0 - upwindWeight) * insideMomentum); }(); - return transportingVelocity * transportedMomentum * scvf.directionSign() * Extrusion::area(scvf) * extrusionFactor_(elemVolVars, scvf); + return transportingVelocity * transportedMomentum * scvf.directionSign() * Extrusion::area(fvGeometry, scvf) * extrusionFactor_(elemVolVars, scvf); } private: diff --git a/dumux/freeflow/navierstokes/momentum/localresidual.hh b/dumux/freeflow/navierstokes/momentum/localresidual.hh index 43924b431fec86cfa57ecac16fff0459dcb0a489..68cb238c3855c4c5af78b36738baeb88bce3ece1 100644 --- a/dumux/freeflow/navierstokes/momentum/localresidual.hh +++ b/dumux/freeflow/navierstokes/momentum/localresidual.hh @@ -240,7 +240,7 @@ public: if (bcTypes.hasNeumann() && bcTypes.isNeumann(scv.dofAxis())) { const auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); - return neumannFluxes[scv.dofAxis()] * Extrusion::area(scvf) * elemVolVars[scv].extrusionFactor(); + return neumannFluxes[scv.dofAxis()] * Extrusion::area(fvGeometry, scvf) * elemVolVars[scv].extrusionFactor(); } } else if (scvf.isLateral()) @@ -273,7 +273,7 @@ public: if (bcTypes.hasNeumann() && bcTypes.isNeumann(scvf.normalAxis())) { const auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); - return neumannFluxes[scv.dofAxis()] * Extrusion::area(scvf) * elemVolVars[scv].extrusionFactor(); + return neumannFluxes[scv.dofAxis()] * Extrusion::area(fvGeometry, scvf) * elemVolVars[scv].extrusionFactor(); } } } @@ -303,7 +303,7 @@ public: if (bcTypes.hasNeumann() && bcTypes.isNeumann(scv.dofAxis())) { const auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); - return neumannFluxes[scv.dofAxis()] * Extrusion::area(scvf) * elemVolVars[scv].extrusionFactor(); + return neumannFluxes[scv.dofAxis()] * Extrusion::area(fvGeometry, scvf) * elemVolVars[scv].extrusionFactor(); } } diff --git a/dumux/freeflow/navierstokes/scalarfluxvariables.hh b/dumux/freeflow/navierstokes/scalarfluxvariables.hh index 3a07b5e9bba8e17d70318e97e84f0b231ca21624..f440ca34bada9971c1e4b54b6ccdc72881949fcf 100644 --- a/dumux/freeflow/navierstokes/scalarfluxvariables.hh +++ b/dumux/freeflow/navierstokes/scalarfluxvariables.hh @@ -69,7 +69,7 @@ public: const ElementFluxVariablesCache& elemFluxVarsCache) { const auto velocity = problem.faceVelocity(element, fvGeometry, scvf); - const Scalar volumeFlux = velocity*scvf.unitOuterNormal()*Extrusion::area(scvf)*extrusionFactor_(elemVolVars, scvf); + const Scalar volumeFlux = velocity*scvf.unitOuterNormal()*Extrusion::area(fvGeometry, scvf)*extrusionFactor_(elemVolVars, scvf); return volumeFlux; } }; @@ -84,7 +84,7 @@ public: { const auto& scvf = this->scvFace(); const auto velocity = this->problem().faceVelocity(this->element(), this->fvGeometry(), scvf); - const Scalar volumeFlux = velocity*scvf.unitOuterNormal()*Extrusion::area(scvf)*extrusionFactor_(this->elemVolVars(), scvf); + const Scalar volumeFlux = velocity*scvf.unitOuterNormal()*Extrusion::area(this->fvGeometry(), scvf)*extrusionFactor_(this->elemVolVars(), scvf); return UpwindScheme::apply(*this, upwindTerm, volumeFlux, 0/*phaseIdx*/); } else diff --git a/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh b/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh index 219bab9561977df3fe42f490092cb69d039d8f94..3873fc1564db00458e6b9884887796667c69d9c9 100644 --- a/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh +++ b/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh @@ -258,7 +258,7 @@ public: const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; const Scalar extrusionFactor = harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor()); - result[0] = elemFaceVars[scvf].velocitySelf() * Extrusion::area(scvf) * extrusionFactor * scvf.directionSign(); + result[0] = elemFaceVars[scvf].velocitySelf() * Extrusion::area(fvGeometry, scvf) * extrusionFactor * scvf.directionSign(); return result; }; diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh index 2ab1a75a9dd0b917a91600ffd3d4fe68fcb6daa3..9f04ee91af8ca6f5890fb0710491a7e3e9d38046 100644 --- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh @@ -112,6 +112,7 @@ public: * \param upwindTerm The uwind term (i.e. the advectively transported quantity) */ template + [[deprecated("Will be removed after release 3.6. Use interface with additional fvGeometry parameter instead.")]] static Scalar advectiveFluxForCellCenter(const Problem& problem, const ElementVolumeVariables& elemVolVars, const ElementFaceVariables& elemFaceVars, @@ -135,6 +136,40 @@ public: return flux * extrusionFactor_(elemVolVars, scvf); } + /*! + * \brief Returns the advective flux over a sub control volume face. + * \param problem The object specifying the problem which ought to be simulated + * \param fvGeometry The element-centered finite volume geometry view + * \param elemVolVars All volume variables for the element + * \param elemFaceVars The face variables + * \param scvf The sub control volume face + * \param upwindTerm The uwind term (i.e. the advectively transported quantity) + */ + template + static Scalar advectiveFluxForCellCenter(const Problem& problem, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const ElementFaceVariables& elemFaceVars, + const SubControlVolumeFace &scvf, + UpwindTerm upwindTerm) + { + const Scalar velocity = elemFaceVars[scvf].velocitySelf(); + const bool insideIsUpstream = scvf.directionSign() == sign(velocity); + static const Scalar upwindWeight = getParamFromGroup(problem.paramGroup(), "Flux.UpwindWeight"); + + const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; + const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; + + const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars; + const auto& downstreamVolVars = insideIsUpstream ? outsideVolVars : insideVolVars; + + const Scalar flux = (upwindWeight * upwindTerm(upstreamVolVars) + + (1.0 - upwindWeight) * upwindTerm(downstreamVolVars)) + * velocity * Extrusion::area(fvGeometry, scvf) * scvf.directionSign(); + + return flux * extrusionFactor_(elemVolVars, scvf); + } + /*! * \brief Computes the flux for the cell center residual (mass balance). * @@ -163,7 +198,7 @@ public: // Call the generic flux function. CellCenterPrimaryVariables result(0.0); - result[Indices::conti0EqIdx - ModelTraits::dim()] = advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTerm); + result[Indices::conti0EqIdx - ModelTraits::dim()] = advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTerm); return result; } @@ -255,7 +290,7 @@ public: // our velocity dof of interest lives on but with adjusted centroid const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); FaceFrontalSubControlVolumeFace frontalFace(scv.center(), scvf.area()); - return frontalFlux * Extrusion::area(frontalFace) * insideVolVars.extrusionFactor(); + return frontalFlux * Extrusion::area(fvGeometry, frontalFace) * insideVolVars.extrusionFactor(); } /*! @@ -334,7 +369,7 @@ public: const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); const auto lateralBoundaryFace = makeStaggeredBoundaryFace(scvf, lateralStaggeredFaceCenter); lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(lateralFace.directionIndex())] - * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(lateralScvf) * lateralFace.directionSign(); + * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(fvGeometry, lateralScvf) * lateralFace.directionSign(); continue; } @@ -350,7 +385,7 @@ public: const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralFace, lateralStaggeredFaceCenter); lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())] - * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(lateralScvf) * lateralFace.directionSign(); + * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(fvGeometry, lateralScvf) * lateralFace.directionSign(); continue; } @@ -374,7 +409,7 @@ public: const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralFace, lateralStaggeredFaceCenter); lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())] - * elemVolVars[lateralFace.insideScvIdx()].extrusionFactor() * Extrusion::area(lateralScvf); + * elemVolVars[lateralFace.insideScvIdx()].extrusionFactor() * Extrusion::area(fvGeometry, lateralScvf); continue; } @@ -431,6 +466,7 @@ public: * // boundary * \endverbatim */ + [[deprecated("Will be removed after release 3.6. Use interface with additional fvGeometry parameter instead.")]] FacePrimaryVariables inflowOutflowBoundaryFlux(const Problem& problem, const Element& element, const SubControlVolumeFace& scvf, @@ -462,6 +498,54 @@ public: return inOrOutflow * scvf.directionSign() * Extrusion::area(scvf) * insideVolVars.extrusionFactor(); } + /*! + * \brief Returns the momentum flux over an inflow or outflow boundary face. + * + * \verbatim + * scvf // + * ---------=======// == and # staggered half-control-volume + * | || #// current scvf + * | || #// # staggered boundary face over which fluxes are calculated + * | || x~~~~> vel.Self + * | || #// x dof position + * scvf | || #// + * --------========// -- element + * scvf // + * // boundary + * \endverbatim + */ + FacePrimaryVariables inflowOutflowBoundaryFlux(const Problem& problem, + const SubControlVolumeFace& scvf, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const ElementFaceVariables& elemFaceVars) const + { + FacePrimaryVariables inOrOutflow(0.0); + const auto& element = fvGeometry.element(); + const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; + + // Advective momentum flux. + if (problem.enableInertiaTerms()) + { + const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf(); + const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; + const auto& upVolVars = (scvf.directionSign() == sign(velocitySelf)) ? + insideVolVars : outsideVolVars; + + inOrOutflow += velocitySelf * velocitySelf * upVolVars.density(); + } + + // Apply a pressure at the boundary. + const Scalar boundaryPressure = normalizePressure + ? (problem.dirichlet(element, scvf)[Indices::pressureIdx] - + problem.initial(scvf)[Indices::pressureIdx]) + : problem.dirichlet(element, scvf)[Indices::pressureIdx]; + inOrOutflow += boundaryPressure; + + // Account for the orientation of the face at the boundary, + return inOrOutflow * scvf.directionSign() * Extrusion::area(fvGeometry, scvf) * insideVolVars.extrusionFactor(); + } + private: /*! @@ -533,7 +617,7 @@ private: StaggeredUpwindHelper upwindHelper(element, fvGeometry, scvf, elemFaceVars, elemVolVars, gridFluxVarsCache.staggeredUpwindMethods()); FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area()); return upwindHelper.computeUpwindLateralMomentum(selfIsUpstream, lateralFace, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes) - * transportingVelocity * lateralFace.directionSign() * Extrusion::area(lateralScvf) * extrusionFactor_(elemVolVars, lateralFace); + * transportingVelocity * lateralFace.directionSign() * Extrusion::area(fvGeometry, lateralScvf) * extrusionFactor_(elemVolVars, lateralFace); } /*! @@ -611,7 +695,7 @@ private: // Account for the area of the staggered lateral face (0.5 of the coinciding scfv). FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area()); - return lateralDiffusiveFlux * Extrusion::area(lateralScvf) * extrusionFactor_(elemVolVars, lateralFace); + return lateralDiffusiveFlux * Extrusion::area(fvGeometry, lateralScvf) * extrusionFactor_(elemVolVars, lateralFace); } /*! @@ -688,7 +772,7 @@ private: FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area()); const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralFace, lateralStaggeredFaceCenter_(scvf, localSubFaceIdx)); lateralFlux += problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())] - * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(lateralScvf); + * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(fvGeometry, lateralScvf); return true; } else diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh index 55d905da705a01bea4ec5a10af241e96a1d5d34f..62adae2bd7bd9ee3c6ab8ae06a7dff9dfb295c15 100644 --- a/dumux/freeflow/navierstokes/staggered/localresidual.hh +++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh @@ -253,7 +253,7 @@ public: for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx) { if (bcTypes.isNeumann(eqIdx + cellCenterOffset)) - result[eqIdx] = neumannFluxes[eqIdx + cellCenterOffset] * extrusionFactor * Extrusion::area(scvf); + result[eqIdx] = neumannFluxes[eqIdx + cellCenterOffset] * extrusionFactor * Extrusion::area(fvGeometry, scvf); } } @@ -325,7 +325,7 @@ public: // add a given Neumann flux for the face on the boundary itself ... const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor(); result = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())] - * extrusionFactor * Extrusion::area(scvf); + * extrusionFactor * Extrusion::area(fvGeometry, scvf); // ... and treat the fluxes of the remaining (frontal and lateral) faces of the staggered control volume result += fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache()); @@ -337,7 +337,7 @@ public: result = fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache()); // incorporate the inflow or outflow contribution - result += fluxVars.inflowOutflowBoundaryFlux(problem, element, scvf, elemVolVars, elemFaceVars); + result += fluxVars.inflowOutflowBoundaryFlux(problem, scvf, fvGeometry, elemVolVars, elemFaceVars); } } return result; @@ -370,7 +370,7 @@ private: if(problem.useWallFunction(element, scvf, eqIdx + cellCenterOffset)) { boundaryFlux[eqIdx] = problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx] - * extrusionFactor * Extrusion::area(scvf); + * extrusionFactor * Extrusion::area(fvGeometry, scvf); } } } diff --git a/dumux/freeflow/nonisothermal/localresidual.hh b/dumux/freeflow/nonisothermal/localresidual.hh index 2a5bc519a9f89bfd75ec4158ab6ca2e4f4ecda0c..d972ae00841c9c15b2bc9323d1caa7b693c7a72f 100644 --- a/dumux/freeflow/nonisothermal/localresidual.hh +++ b/dumux/freeflow/nonisothermal/localresidual.hh @@ -104,6 +104,7 @@ public: auto upwindTerm = [](const auto& volVars) { return volVars.density() * volVars.enthalpy(); }; flux[localEnergyBalanceIdx] += FluxVariables::advectiveFluxForCellCenter(problem, + fvGeometry, elemVolVars, elemFaceVars, scvf, diff --git a/dumux/freeflow/rans/oneeq/staggered/fluxvariables.hh b/dumux/freeflow/rans/oneeq/staggered/fluxvariables.hh index 8b18e0d18aa253f96bf07105741dc886a0bd6ee0..0c4289a52166ed61bf6768b39885b2eb5071e663 100644 --- a/dumux/freeflow/rans/oneeq/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/oneeq/staggered/fluxvariables.hh @@ -100,7 +100,7 @@ public: }; flux[viscosityTildeEqIdx] - = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermK); + = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTermK); // calculate diffusive flux const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); @@ -142,7 +142,7 @@ public: flux[viscosityTildeEqIdx] += coeff / distance * (insideVolVars.viscosityTilde() - outsideVolVars.viscosityTilde()) - * Extrusion::area(scvf); + * Extrusion::area(fvGeometry, scvf); } return flux; } diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh index a7232f9e056891ed00c6fa3a7aada4c0b73db083..2a297d3a4b793185d3199a3cb458774ee4e7a563 100644 --- a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh @@ -106,9 +106,9 @@ public: }; flux[turbulentKineticEnergyEqIdx] - = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermK); + = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTermK); flux[dissipationEqIdx ] - = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermEpsilon); + = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTermEpsilon); // calculate diffusive flux const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); @@ -159,7 +159,7 @@ public: flux[turbulentKineticEnergyEqIdx] += coeff_k / distance * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) - * Extrusion::area(scvf); + * Extrusion::area(fvGeometry, scvf); } } @@ -169,7 +169,7 @@ public: flux[dissipationEqIdx] += coeff_e / distance * (insideVolVars.dissipation() - outsideVolVars.dissipation()) - * Extrusion::area(scvf); + * Extrusion::area(fvGeometry, scvf); } return flux; } @@ -190,7 +190,7 @@ public: return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) + 2.0 / ModelTraits::dim() * insideVolVars.density() * insideVolVars.turbulentKineticEnergy() - * Extrusion::area(scvf) * scvf.directionSign() * insideVolVars.extrusionFactor(); + * Extrusion::area(fvGeometry, scvf) * scvf.directionSign() * insideVolVars.extrusionFactor(); } }; diff --git a/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh index e8b86e24b2370f756985bb332d8335e278083bd2..e347aa55d7b0837e0fc0a8a82c41b66f2d38fe54 100644 --- a/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh @@ -105,9 +105,9 @@ public: }; flux[turbulentKineticEnergyEqIdx] - = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermK); + = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTermK); flux[dissipationEqIdx] - = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermOmega); + = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTermOmega); // calculate diffusive flux const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); @@ -159,7 +159,7 @@ public: flux[turbulentKineticEnergyEqIdx] += coeff_k / distance * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) - * Extrusion::area(scvf); + * Extrusion::area(fvGeometry, scvf); } if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::dissipationEqIdx) || bcTypes.isSymmetry()))) @@ -167,7 +167,7 @@ public: flux[dissipationEqIdx] += coeff_w / distance * (insideVolVars.dissipation() - outsideVolVars.dissipation()) - * Extrusion::area(scvf); + * Extrusion::area(fvGeometry, scvf); } return flux; } @@ -188,7 +188,7 @@ public: return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) + 2.0 / ModelTraits::dim() * insideVolVars.density() * insideVolVars.turbulentKineticEnergy() - * Extrusion::area(scvf) * scvf.directionSign() * insideVolVars.extrusionFactor(); + * Extrusion::area(fvGeometry, scvf) * scvf.directionSign() * insideVolVars.extrusionFactor(); } }; diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh index b6aa6f908622b3ad96f759c8b0a7f85c48aef9e0..874dd54aef364dfea36b5d6f8d9565de48fabc53 100644 --- a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh @@ -105,9 +105,9 @@ public: }; flux[turbulentKineticEnergyEqIdx] - = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermK); + = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTermK); flux[dissipationEqIdx] - = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermEpsilon); + = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTermEpsilon); // calculate diffusive flux const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); @@ -159,7 +159,7 @@ public: flux[turbulentKineticEnergyEqIdx] += coeff_k / distance * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) - * Extrusion::area(scvf); + * Extrusion::area(fvGeometry, scvf); } if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::dissipationEqIdx) || bcTypes.isSymmetry()))) @@ -167,7 +167,7 @@ public: flux[dissipationEqIdx] += coeff_e / distance * (insideVolVars.dissipationTilde() - outsideVolVars.dissipationTilde()) - * Extrusion::area(scvf); + * Extrusion::area(fvGeometry, scvf); } return flux; } @@ -188,7 +188,7 @@ public: return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) + 2.0 / ModelTraits::dim() * insideVolVars.density() * insideVolVars.turbulentKineticEnergy() - * Extrusion::area(scvf) * scvf.directionSign() * insideVolVars.extrusionFactor(); + * Extrusion::area(fvGeometry, scvf) * scvf.directionSign() * insideVolVars.extrusionFactor(); } }; diff --git a/dumux/freeflow/rans/twoeq/sst/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/sst/staggered/fluxvariables.hh index d6f62f6146cb5949c3fb03dda964aa72a3973b2d..8384e2b485a415f35fe11ab799c352de9128b52c 100644 --- a/dumux/freeflow/rans/twoeq/sst/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/twoeq/sst/staggered/fluxvariables.hh @@ -103,9 +103,9 @@ public: { return volVars.dissipation() * volVars.density(); }; flux[turbulentKineticEnergyEqIdx] - = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermK); + = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTermK); flux[dissipationEqIdx] - = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermOmega); + = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTermOmega); // calculate diffusive flux const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); @@ -174,7 +174,7 @@ public: flux[turbulentKineticEnergyEqIdx] += coeff_k / distance * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) - * Extrusion::area(scvf); + * Extrusion::area(fvGeometry, scvf); } if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::dissipationEqIdx) || bcTypes.isSymmetry()))) @@ -182,7 +182,7 @@ public: flux[dissipationEqIdx] += coeff_w / distance * (insideVolVars.dissipation() - outsideVolVars.dissipation()) - * Extrusion::area(scvf); + * Extrusion::area(fvGeometry, scvf); } return flux; } @@ -203,7 +203,7 @@ public: return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) + 2.0 / ModelTraits::dim() * insideVolVars.density() * insideVolVars.turbulentKineticEnergy() - * Extrusion::area(scvf) * scvf.directionSign() * insideVolVars.extrusionFactor(); + * Extrusion::area(fvGeometry, scvf) * scvf.directionSign() * insideVolVars.extrusionFactor(); } }; diff --git a/dumux/geometry/volume.hh b/dumux/geometry/volume.hh index 9650d2ca8ba657c9a33be7e7546718453e875962..ea2888714fe29bed9341b12188c1f2b5fe0cfc5f 100644 --- a/dumux/geometry/volume.hh +++ b/dumux/geometry/volume.hh @@ -170,13 +170,30 @@ auto convexPolytopeVolume(const Geometry& geo) template auto volume(const Geometry& geo, unsigned int integrationOrder = 4) { - double volume = 0.0; - const auto rule = Dune::QuadratureRules::rule(geo.type(), integrationOrder); + using ctype = typename Geometry::ctype; + ctype volume = 0.0; + const auto rule = Dune::QuadratureRules::rule(geo.type(), integrationOrder); for (const auto& qp : rule) volume += geo.integrationElement(qp.position())*qp.weight(); return volume; } +/*! + * \ingroup Geometry + * \brief The volume of a given geometry with an extrusion/transformation policy + * \note depending on the transformation this might not be an accurate quadrature rule anymore + */ +template +auto volume(const Geometry& geo, Transformation transformation, unsigned int integrationOrder = 4) +{ + using ctype = typename Geometry::ctype; + ctype volume = 0.0; + const auto rule = Dune::QuadratureRules::rule(geo.type(), integrationOrder); + for (const auto& qp : rule) + volume += transformation.integrationElement(geo, qp.position())*qp.weight(); + return volume; +} + } // end namespace Dumux #endif diff --git a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh index 861054043920a749ab66844da3deb3c650ed89aa..2fafb0196089f21e88239a9952be2c1d6c3677b9 100644 --- a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh +++ b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh @@ -178,11 +178,11 @@ public: const auto& insideVolVars = elemVolVars[insideScv]; using Extrusion = typename GridGeometry::Extrusion; - const auto ti = computeTpfaTransmissibility(scvf, insideScv, insideVolVars.permeability(), insideVolVars.extrusionFactor()); - const auto tj = computeTpfaTransmissibility(flipScvf, outsideScv, outsideVolVars.permeability(), outsideVolVars.extrusionFactor()); + const auto ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, insideVolVars.permeability(), insideVolVars.extrusionFactor()); + const auto tj = computeTpfaTransmissibility(fvGeometryOutside, flipScvf, outsideScv, outsideVolVars.permeability(), outsideVolVars.extrusionFactor()); Scalar tij = 0.0; if (ti*tj > 0.0) - tij = Extrusion::area(scvf)*(ti * tj)/(ti + tj); + tij = Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj); if (enableGravity) { @@ -198,7 +198,7 @@ public: const auto pInside = insideVolVars.pressure(0); const auto pOutside = outsideVolVars.pressure(0); - flux = tij*(pInside - pOutside) + rho*Extrusion::area(scvf)*alpha_inside - rho*tij/tj*(alpha_inside - alpha_outside); + flux = tij*(pInside - pOutside) + rho*Extrusion::area(fvGeometry, scvf)*alpha_inside - rho*tij/tj*(alpha_inside - alpha_outside); } else @@ -222,7 +222,7 @@ public: + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); // make this a area-specific flux - flux /= Extrusion::area(scvf)*insideVolVars.extrusionFactor(); + flux /= Extrusion::area(fvGeometry, scvf)*insideVolVars.extrusionFactor(); return flux; } diff --git a/dumux/multidomain/boundary/freeflowporenetwork/couplingmapper.hh b/dumux/multidomain/boundary/freeflowporenetwork/couplingmapper.hh index c8d5e17e23e655482a0e8254f09992bb824856d1..6b10c18b71293b5d76b288b39168cae21c6bdbca 100644 --- a/dumux/multidomain/boundary/freeflowporenetwork/couplingmapper.hh +++ b/dumux/multidomain/boundary/freeflowporenetwork/couplingmapper.hh @@ -325,7 +325,7 @@ private: using std::abs; for (const auto& scv : scvs(fvGeometry)) { - const Scalar eps = diameter(scv.geometry())*1e-6; // TODO + const Scalar eps = diameter(fvGeometry.geometry(scv))*1e-6; // TODO assert(eps < couplingPoreRadius); if (scv.dofAxis() == couplingInterfaceDirectionIdx) // the free flow dofs that lie within the coupling interface diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh index 6603b582688e2a248e125e9f9daf23f3ebc7d094..74e8145786e9c2cb4cc45288013498e942bd7ce5 100644 --- a/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh +++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh @@ -341,7 +341,7 @@ protected: const auto alpha = vtmv(scvf.unitOuterNormal(), K, gravity); const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); - const auto ti = computeTpfaTransmissibility(scvf, insideScv, K, 1.0); + const auto ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, K, 1.0); return (-1/couplingPhaseMobility * (scvf.unitOuterNormal() * couplingPhaseVelocity) + couplingPhaseDensity * alpha)/ti + couplingPhaseCellCenterPressure; diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper.hh index 5bca370ef4b658f47a6ba03f4bd6d82cb4fd4959..1e64decf0270798a90d683392298e7d8e2788698 100644 --- a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper.hh +++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper.hh @@ -130,7 +130,7 @@ public: // get elements intersecting with the scvf center // for robustness add epsilon in unit outer normal direction - const auto eps = (pmScvf.ipGlobal() - pmScvf.geometry().corner(0)).two_norm()*1e-8; + const auto eps = (pmScvf.ipGlobal() - pmFvGeometry.geometry(pmScvf).corner(0)).two_norm()*1e-8; auto globalPos = pmScvf.ipGlobal(); globalPos.axpy(eps, pmScvf.unitOuterNormal()); const auto indices = intersectingEntities(globalPos, freeFlowMomentumGG.boundingBoxTree()); @@ -181,7 +181,7 @@ public: else { // for robustness add epsilon in unit outer normal direction - const auto otherScvfEps = (otherFfScvf.ipGlobal() - otherFfScvf.geometry().corner(0)).two_norm()*1e-8; + const auto otherScvfEps = (otherFfScvf.ipGlobal() - ffFvGeometry.geometry(otherFfScvf).corner(0)).two_norm()*1e-8; auto otherScvfGlobalPos = otherFfScvf.center(); otherScvfGlobalPos.axpy(otherScvfEps, otherFfScvf.unitOuterNormal()); if (!intersectingEntities(otherScvfGlobalPos, porousMediumGG.boundingBoxTree()).empty()) diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh index c768126d90120e71796af8361e7e91321a91e87b..14221ca9a256bc6cff040b0dd3ae4ffdaef486a4 100644 --- a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh @@ -489,7 +489,7 @@ protected: const auto alpha = vtmv(scvf.unitOuterNormal(), K, couplingManager_.problem(darcyIdx).spatialParams().gravity(scvf.center())); const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); - const auto ti = computeTpfaTransmissibility(scvf, insideScv, K, 1.0); + const auto ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, K, 1.0); return (-1/couplingPhaseMobility * (scvf.unitOuterNormal() * couplingPhaseVelocity) + couplingPhaseDensity * alpha)/ti + couplingPhaseCellCenterPressure; diff --git a/dumux/multidomain/embedded/couplingmanagerbase.hh b/dumux/multidomain/embedded/couplingmanagerbase.hh index ff47938755f74158bbc6774ec15d5da924dd70a6..b691cbac07db165588d79ad93c4039d46eb875f2 100644 --- a/dumux/multidomain/embedded/couplingmanagerbase.hh +++ b/dumux/multidomain/embedded/couplingmanagerbase.hh @@ -218,7 +218,7 @@ public: { auto couplingSource = this->problem(domainI).scvPointSources(element, fvGeometry, curElemVolVars, scv); couplingSource += this->problem(domainI).source(element, fvGeometry, curElemVolVars, scv); - couplingSource *= -GridGeometry::Extrusion::volume(scv)*curElemVolVars[scv].extrusionFactor(); + couplingSource *= -GridGeometry::Extrusion::volume(fvGeometry, scv)*curElemVolVars[scv].extrusionFactor(); residual[scv.indexInElement()] = couplingSource; } return residual; diff --git a/dumux/multidomain/facet/box/darcyslaw.hh b/dumux/multidomain/facet/box/darcyslaw.hh index 45b460a338a85f4b5e80e282834806a83a1d4e6e..d404520cf6aba858116d816fa1225c245ba343f9 100644 --- a/dumux/multidomain/facet/box/darcyslaw.hh +++ b/dumux/multidomain/facet/box/darcyslaw.hh @@ -122,7 +122,7 @@ public: gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center())); // apply facet permeability and return the flux - return -1.0*Extrusion::area(scvf) + return -1.0*Extrusion::area(fvGeometry, scvf) *insideVolVars.extrusionFactor() *vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), gradP); } @@ -154,7 +154,7 @@ public: gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center())); // apply matrix permeability and return the flux - return -1.0*Extrusion::area(scvf) + return -1.0*Extrusion::area(fvGeometry, scvf) *insideVolVars.extrusionFactor() *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), gradP); } diff --git a/dumux/multidomain/facet/box/fickslaw.hh b/dumux/multidomain/facet/box/fickslaw.hh index 075aa5f17520f2a376f3c15a8ab073f600c57121..118107a84bf0295ccc8438cc7707ace83ce4eadc 100644 --- a/dumux/multidomain/facet/box/fickslaw.hh +++ b/dumux/multidomain/facet/box/fickslaw.hh @@ -132,7 +132,7 @@ public: auto gradX = preGradX; gradX *= (massOrMoleFraction(facetVolVars, referenceSystem, phaseIdx, compIdx) - x); - componentFlux[compIdx] = -1.0*rho*Extrusion::area(scvf) + componentFlux[compIdx] = -1.0*rho*Extrusion::area(fvGeometry, scvf) *insideVolVars.extrusionFactor() *vtmv(scvf.unitOuterNormal(), facetVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx), @@ -170,7 +170,7 @@ public: gradX.axpy(xFractions[scv.localDofIndex()], fluxVarCache.gradN(scv.indexInElement())); // apply matrix diffusion coefficient and return the flux - componentFlux[compIdx] = -1.0*rho*Extrusion::area(scvf) + componentFlux[compIdx] = -1.0*rho*Extrusion::area(fvGeometry, scvf) *insideVolVars.extrusionFactor() *vtmv(scvf.unitOuterNormal(), insideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx), diff --git a/dumux/multidomain/facet/box/fourierslaw.hh b/dumux/multidomain/facet/box/fourierslaw.hh index a6363228c35c4ab4448f953534fe9227a5e43ccc..5bc2432b01de7364c370e89962d3540d2a59711b 100644 --- a/dumux/multidomain/facet/box/fourierslaw.hh +++ b/dumux/multidomain/facet/box/fourierslaw.hh @@ -117,7 +117,7 @@ public: gradT /= gradT.two_norm2(); gradT *= (facetVolVars.temperature() - T); - return -1.0*Extrusion::area(scvf) + return -1.0*Extrusion::area(fvGeometry, scvf) *insideVolVars.extrusionFactor() *vtmv(scvf.unitOuterNormal(), facetVolVars.effectiveThermalConductivity(), @@ -144,7 +144,7 @@ public: gradT.axpy(temperatures[scv.localDofIndex()], fluxVarCache.gradN(scv.indexInElement())); // apply matrix thermal conductivity and return the flux - return -1.0*Extrusion::area(scvf) + return -1.0*Extrusion::area(fvGeometry, scvf) *insideVolVars.extrusionFactor() *vtmv(scvf.unitOuterNormal(), insideVolVars.effectiveThermalConductivity(), diff --git a/dumux/multidomain/facet/box/localresidual.hh b/dumux/multidomain/facet/box/localresidual.hh index 7d1d2568e97217bb36d55a1290216ab2c05c7c99..cfba0bbcfc74c74a98d303fbc1a588638707d504 100644 --- a/dumux/multidomain/facet/box/localresidual.hh +++ b/dumux/multidomain/facet/box/localresidual.hh @@ -116,7 +116,7 @@ public: auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); // multiply neumann fluxes with the area and the extrusion factor - neumannFluxes *= Extrusion::area(scvf)*elemVolVars[scv].extrusionFactor(); + neumannFluxes *= Extrusion::area(fvGeometry, scvf)*elemVolVars[scv].extrusionFactor(); // only add fluxes to equations for which Neumann is set for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx) diff --git a/dumux/multidomain/facet/box/subcontrolvolumeface.hh b/dumux/multidomain/facet/box/subcontrolvolumeface.hh index 207de2a27bc398e4fe2649acc6263a86900f9184..c3b14ec4bddb853a2f3dccdd26f0cf7ee539dfda 100644 --- a/dumux/multidomain/facet/box/subcontrolvolumeface.hh +++ b/dumux/multidomain/facet/box/subcontrolvolumeface.hh @@ -35,6 +35,7 @@ #include #include #include +#include namespace Dumux { @@ -78,7 +79,10 @@ public: : corners_(geometryHelper.getScvfCorners(scvfIndex)) , center_(0.0) , unitOuterNormal_(geometryHelper.normal(corners_, scvIndices)) - , area_(geometryHelper.scvfArea(corners_)) + , area_(Dumux::convexPolytopeVolume( + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; }) + ) , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) , facetIndex_(/*undefined*/) @@ -105,7 +109,10 @@ public: : corners_(geometryHelper.getBoundaryScvfCorners(intersection.indexInInside(), indexInIntersection)) , center_(0.0) , unitOuterNormal_(intersection.centerUnitOuterNormal()) - , area_(geometryHelper.scvfArea(corners_)) + , area_(Dumux::convexPolytopeVolume( + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; }) + ) , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) , facetIndex_(intersection.indexInInside()) diff --git a/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh b/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh index d28b08cb6b30124b88a35858b358d32d3cbd46ec..61640013b26de3f1a5edfe1d77b4754f9b33b2a7 100644 --- a/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh +++ b/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh @@ -229,7 +229,7 @@ public: rho /= numOutsideFaces + 1; deltaG[localDofIdx] -= curXiFactor*alpha_inside; - deltaG[localDofIdx] *= rho*Extrusion::area(curGlobalScvf); + deltaG[localDofIdx] *= rho*Extrusion::area(this->fvGeometry(), curGlobalScvf); } // use average density between facet and cell density on interior Dirichlet boundaries else if (curIsInteriorBoundary) @@ -243,13 +243,13 @@ public: rho = getRho(this->elemVolVars()[curGlobalScvf.outsideScvIdx()]); // add "inside" & "outside" alphas to gravity containers - g[faceIdx] = alpha_inside*rho*Extrusion::area(curGlobalScvf); + g[faceIdx] = alpha_inside*rho*Extrusion::area(this->fvGeometry(), curGlobalScvf); if (isSurfaceGrid && !curIsInteriorBoundary) { unsigned int i = 0; for (const auto& alpha : alpha_outside) - outsideG[faceIdx][i++] = alpha*rho*Extrusion::area(curGlobalScvf); + outsideG[faceIdx][i++] = alpha*rho*Extrusion::area(this->fvGeometry(), curGlobalScvf); } } @@ -340,7 +340,7 @@ private: // the omega factors of the "positive" sub volume Helper::resizeVector(wijk[faceIdx], /*no outside scvs present*/1); - wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); + wijk[faceIdx][0] = computeMpfaTransmissibility(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); for (LocalIndexType localDir = 0; localDir < dim; localDir++) @@ -379,7 +379,7 @@ private: // the omega factors of the "positive" sub volume Helper::resizeVector(wijk[faceIdx], neighborScvIndices.size()); - wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); + wijk[faceIdx][0] = computeMpfaTransmissibility(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); using std::abs; bool isZeroWij = false; @@ -434,7 +434,7 @@ private: // On surface grids we use the square root of the extrusion factor as approximation of the aperture using std::sqrt; - const auto wFacet = 2.0*Extrusion::area(curGlobalScvf)*posVolVars.extrusionFactor() + const auto wFacet = 2.0*Extrusion::area(this->fvGeometry(), curGlobalScvf)*posVolVars.extrusionFactor() *vtmv(curGlobalScvf.unitOuterNormal(), facetTensor, curGlobalScvf.unitOuterNormal()) / (dim < dimWorld ? sqrt(facetVolVars.extrusionFactor()) : facetVolVars.extrusionFactor()); @@ -467,7 +467,7 @@ private: // On surface grids, use outside face for "negative" transmissibility calculation const auto& scvf = dim < dimWorld ? this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside) : curGlobalScvf; - wijk[faceIdx][idxOnScvf] = computeMpfaTransmissibility(negLocalScv, scvf, negTensor, negVolVars.extrusionFactor()); + wijk[faceIdx][idxOnScvf] = computeMpfaTransmissibility(this->fvGeometry(), negLocalScv, scvf, negTensor, negVolVars.extrusionFactor()); // flip sign on surface grids (since we used the "outside" normal) if (dim < dimWorld) diff --git a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh index 05c4ec5d7abb274ae2604eda9af6ccbfd265baec..bd7ca6dbfa4afa9afa73d0e6d04a6319796f2085 100644 --- a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh +++ b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh @@ -193,7 +193,7 @@ class CCTpfaFacetCouplingDarcysLawImpllocalResidual().computeSource(problem(), element, this->fvGeometry(), elemVolVars, scv); - source *= -Extrusion::volume(scv)*curVolVars.extrusionFactor(); + source *= -Extrusion::volume(this->fvGeometry(), scv)*curVolVars.extrusionFactor(); residual[scv.localDofIndex()] = std::move(source); } diff --git a/dumux/multidomain/subdomaincclocalassembler.hh b/dumux/multidomain/subdomaincclocalassembler.hh index fb79b537fc93c260ed74df7af59ebb1d527ac77e..5cbb38c64b0d83b948cacb2c499bd562c5264985 100644 --- a/dumux/multidomain/subdomaincclocalassembler.hh +++ b/dumux/multidomain/subdomaincclocalassembler.hh @@ -190,7 +190,7 @@ public: { const auto& curVolVars = elemVolVars[scv]; auto source = this->localResidual().computeSource(problem(), element, this->fvGeometry(), elemVolVars, scv); - source *= -Extrusion::volume(scv)*curVolVars.extrusionFactor(); + source *= -Extrusion::volume(this->fvGeometry(), scv)*curVolVars.extrusionFactor(); residual[scv.indexInElement()] = std::move(source); } diff --git a/dumux/porenetwork/2p/elementfluxvariablescache.hh b/dumux/porenetwork/2p/elementfluxvariablescache.hh index e111fc1eb95b2c54764268ab61c5025b5bfa4c8d..92bd6b8bbda6e47c0bcbe9fc78a075fb3ea23790 100644 --- a/dumux/porenetwork/2p/elementfluxvariablescache.hh +++ b/dumux/porenetwork/2p/elementfluxvariablescache.hh @@ -28,7 +28,7 @@ #include #include -#include +#include namespace Dumux::PoreNetwork { @@ -44,9 +44,9 @@ class PNMTwoPElementFluxVariablesCache; * \brief The flux variables caches for an element with caching enabled */ template -class PNMTwoPElementFluxVariablesCache : public BoxElementFluxVariablesCache +class PNMTwoPElementFluxVariablesCache : public CVFEElementFluxVariablesCache { - using ParentType = BoxElementFluxVariablesCache; + using ParentType = CVFEElementFluxVariablesCache; public: using ParentType::ParentType; }; diff --git a/dumux/porenetwork/2p/gridfluxvariablescache.hh b/dumux/porenetwork/2p/gridfluxvariablescache.hh index 6b05046147a06a0e4b20da6b8c99f84475c67977..1dcb53ac5cfb7e46d8e4868c7d93e8ba83fd62d5 100644 --- a/dumux/porenetwork/2p/gridfluxvariablescache.hh +++ b/dumux/porenetwork/2p/gridfluxvariablescache.hh @@ -26,7 +26,7 @@ #include #include -#include +#include #include "elementfluxvariablescache.hh" #include "invasionstate.hh" diff --git a/dumux/porousmediumflow/2p/griddatatransfer.hh b/dumux/porousmediumflow/2p/griddatatransfer.hh index 38ec15039d28bed76154fd81bee6111a04ab1b9d..2b35f708afab6cb1e5ae243d8daf29386e009cb8 100644 --- a/dumux/porousmediumflow/2p/griddatatransfer.hh +++ b/dumux/porousmediumflow/2p/griddatatransfer.hh @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -163,7 +164,7 @@ public: VolumeVariables volVars; volVars.update(adaptedValues.u, *problem_, element, scv); - const auto poreVolume = Extrusion::volume(scv)*volVars.porosity(); + const auto poreVolume = Extrusion::volume(fvGeometry, scv)*volVars.porosity(); adaptedValues.associatedMass[phase1Idx] += poreVolume * volVars.density(phase1Idx) * volVars.saturation(phase1Idx); adaptedValues.associatedMass[phase0Idx] += poreVolume * volVars.density(phase0Idx) * volVars.saturation(phase0Idx); } @@ -240,7 +241,7 @@ public: if (!isBox) elemSol[0] /= adaptedValues.count; - const auto elementVolume = Extrusion::volume(element.geometry()); + const auto elementVolume = volume(element.geometry(), Extrusion{}); for (const auto& scv : scvs(fvGeometry)) { VolumeVariables volVars; @@ -272,7 +273,7 @@ public: // For the box scheme, add mass & mass coefficient to container (saturations are recalculated at the end) else { - const auto scvVolume = Extrusion::volume(scv); + const auto scvVolume = Extrusion::volume(fvGeometry, scv); if (formulation == p0s1) { massCoeff[dofIdxGlobal] += scvVolume * volVars.density(phase1Idx) * volVars.porosity(); @@ -324,11 +325,11 @@ public: // overwrite the saturation by a mass conservative one here Scalar massCoeffSon = 0.0; if (formulation == p0s1) - massCoeffSon = Extrusion::volume(scv) * volVars.density(phase1Idx) * volVars.porosity(); + massCoeffSon = Extrusion::volume(fvGeometry, scv) * volVars.density(phase1Idx) * volVars.porosity(); else if (formulation == p1s0) - massCoeffSon = Extrusion::volume(scv) * volVars.density(phase0Idx) * volVars.porosity(); + massCoeffSon = Extrusion::volume(fvGeometry, scv) * volVars.density(phase0Idx) * volVars.porosity(); sol_[scv.dofIndex()][saturationIdx] = - ( Extrusion::volume(scv)/Extrusion::volume(fatherElement.geometry())*massFather )/massCoeffSon; + ( Extrusion::volume(fvGeometry, scv)/volume(fatherElement.geometry(), Extrusion{})*massFather )/massCoeffSon; } } else @@ -347,14 +348,14 @@ public: scv.dofPosition()); // compute mass & mass coefficients for the scvs (saturations are recalculated at the end) - const auto fatherElementVolume = Extrusion::volume(fatherGeometry); + const auto fatherElementVolume = volume(fatherGeometry, Extrusion{}); for (const auto& scv : scvs(fvGeometry)) { VolumeVariables volVars; volVars.update(elemSolSon, *problem_, element, scv); const auto dofIdxGlobal = scv.dofIndex(); - const auto scvVolume = Extrusion::volume(scv); + const auto scvVolume = Extrusion::volume(fvGeometry, scv); if (formulation == p0s1) { massCoeff[dofIdxGlobal] += scvVolume * volVars.density(phase1Idx) * volVars.porosity(); diff --git a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh index 659f70db00b9683a9a079163604771b5efbaf0f9..922fed8670a3866b155adbb94bd9b12ef59e6f1c 100644 --- a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh +++ b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh @@ -105,7 +105,7 @@ public: static_assert(ModelTraits::priVarFormulation() == TwoPFormulation::p0s1, "2p/incompressiblelocalresidual.hh: Analytic differentiation has to be checked for p1-s0 formulation!"); - const auto poreVolume = Extrusion::volume(scv)*curVolVars.porosity(); + const auto poreVolume = Extrusion::volume(fvGeometry, scv)*curVolVars.porosity(); const auto dt = this->timeLoop().timeStepSize(); // partial derivative of the phase storage terms w.r.t. S_n diff --git a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh index cb880d192caef700bd8d968ae1f1ccb87b756f72..fb91e85fc4cb935795febc06a8355784b9ea514e 100644 --- a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh +++ b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh @@ -32,6 +32,7 @@ #include #include #include +#include namespace Dumux { @@ -101,7 +102,10 @@ public: : isFractureScv_(false) , corners_(geometryHelper.getScvCorners(scvIdx)) , center_(0.0) - , volume_(geometryHelper.scvVolume(corners_)) + , volume_(Dumux::convexPolytopeVolume( + Dune::GeometryTypes::cube(T::dim), + [&](unsigned int i){ return corners_[i]; }) + ) , elementIndex_(elementIndex) , vIdxLocal_(scvIdx) , elemLocalScvIdx_(scvIdx) @@ -151,7 +155,10 @@ public: corners_[i] = corners[i]; // compute volume and scv center - volume_ = geometryHelper.scvfArea(corners); + volume_ = Dumux::convexPolytopeVolume( + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; } + ); for (const auto& corner : corners_) center_ += corner; center_ /= corners_.size(); diff --git a/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh index 6c7f0df124d8da93b2c0d81dab2b8e091c54750d..9a1a24a9fa6fbde8bc1a82fbb1ab86d683fbd0d8 100644 --- a/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh +++ b/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh @@ -33,6 +33,7 @@ #include #include #include +#include namespace Dumux { @@ -103,7 +104,10 @@ public: : corners_(geometryHelper.getScvfCorners(scvfIndex)) , center_(0.0) , unitOuterNormal_(geometryHelper.normal(corners_, scvIndices)) - , area_(geometryHelper.scvfArea(corners_)) + , area_(Dumux::convexPolytopeVolume( + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; }) + ) , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) , boundary_(false) @@ -127,7 +131,10 @@ public: : corners_(geometryHelper.getBoundaryScvfCorners(intersection.indexInInside(), indexInIntersection)) , center_(0.0) , unitOuterNormal_(intersection.centerUnitOuterNormal()) - , area_(geometryHelper.scvfArea(corners_)) + , area_(Dumux::convexPolytopeVolume( + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; }) + ) , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) , boundary_(true) diff --git a/dumux/porousmediumflow/fluxvariablescachefiller.hh b/dumux/porousmediumflow/fluxvariablescachefiller.hh index 1ec42a8d138fce83c1da234ec120830f063ebe87..b1a3d10a0215c925fbc509b4a11fd70d761cb0cb 100644 --- a/dumux/porousmediumflow/fluxvariablescachefiller.hh +++ b/dumux/porousmediumflow/fluxvariablescachefiller.hh @@ -637,8 +637,8 @@ private: const auto& scv = fvGeometry().scv(iv.localScv(0).gridScvIndex()); const auto& scvf = fvGeometry().scvf(iv.localScvf(0).gridScvfIndex()); const auto& vv = elemVolVars()[scv]; - const auto eps = Extrusion::area(scvf)*computeTpfaTransmissibility( - scvf, scv, zeroD, vv.extrusionFactor() + const auto eps = Extrusion::area(fvGeometry(), scvf)*computeTpfaTransmissibility( + fvGeometry(), scvf, scv, zeroD, vv.extrusionFactor() ); localAssembler.assembleMatrices(handle.diffusionHandle(), iv, getD, eps); diff --git a/dumux/porousmediumflow/richards/localresidual.hh b/dumux/porousmediumflow/richards/localresidual.hh index 218c06a5bea39c648a4b773482b374c0e26dcf3c..a3fc5ed588bcaf9a13332cf72e2df7a2c4dbd30d 100644 --- a/dumux/porousmediumflow/richards/localresidual.hh +++ b/dumux/porousmediumflow/richards/localresidual.hh @@ -209,7 +209,7 @@ public: static_assert(!FluidSystem::isCompressible(0), "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!"); - const auto poreVolume = Extrusion::volume(scv)*curVolVars.porosity()*curVolVars.extrusionFactor(); + const auto poreVolume = Extrusion::volume(fvGeometry, scv)*curVolVars.porosity()*curVolVars.extrusionFactor(); static const auto rho = curVolVars.density(0); // partial derivative of storage term w.r.t. p_w diff --git a/dumux/porousmediumflow/tracer/localresidual.hh b/dumux/porousmediumflow/tracer/localresidual.hh index 754d5b568c007689506ba3849a9bcef654eb2336..770674990e503633e78d4fe5ea9c032447b0108d 100644 --- a/dumux/porousmediumflow/tracer/localresidual.hh +++ b/dumux/porousmediumflow/tracer/localresidual.hh @@ -218,7 +218,7 @@ public: const auto porosity = curVolVars.porosity(); const auto rho = useMoles ? curVolVars.molarDensity() : curVolVars.density(); - const auto d_storage = Extrusion::volume(scv)*porosity*rho*saturation/this->timeLoop().timeStepSize(); + const auto d_storage = Extrusion::volume(fvGeometry, scv)*porosity*rho*saturation/this->timeLoop().timeStepSize(); for (int compIdx = 0; compIdx < numComponents; ++compIdx) partialDerivatives[compIdx][compIdx] += d_storage; diff --git a/dumux/porousmediumflow/velocity.hh b/dumux/porousmediumflow/velocity.hh index 8f0af5f7022e9bd4ee6a740bdabcc807b3054443..6f89955df46854dc6fa8d6ad61556d394bc7f733 100644 --- a/dumux/porousmediumflow/velocity.hh +++ b/dumux/porousmediumflow/velocity.hh @@ -163,7 +163,7 @@ public: Scalar localArea = scvfReferenceArea_(geomType, scvf.index()); Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea; const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; - flux /= insideVolVars.extrusionFactor() * Extrusion::area(scvf) / scvf.area(); + flux /= insideVolVars.extrusionFactor() * Extrusion::area(fvGeometry, scvf) / scvf.area(); tmpVelocity *= flux; const int eIdxGlobal = gridGeometry_.elementMapper().index(element); @@ -208,7 +208,7 @@ public: Scalar localArea = scvfReferenceArea_(geomType, scvf.index()); Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea; const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; - flux /= insideVolVars.extrusionFactor() * Extrusion::area(scvf) / scvf.area(); + flux /= insideVolVars.extrusionFactor() * Extrusion::area(fvGeometry, scvf) / scvf.area(); // transform the volume flux into a velocity vector Velocity tmpVelocity = localNormal; diff --git a/examples/liddrivencavity/doc/problem.md b/examples/liddrivencavity/doc/problem.md index e70cd12f88b5f66a8bb214caa2677cdb83befd38..e2d607cd1e93a7df35d066ba7aca60853f87dc0f 100644 --- a/examples/liddrivencavity/doc/problem.md +++ b/examples/liddrivencavity/doc/problem.md @@ -498,7 +498,8 @@ This is done for both the momentum and mass grid geometries auto massGridGeometry = std::make_shared(leafGridView); ``` -We introduce the multidomain coupling manager, which will coupled the mass and the momentum problems +We introduce the multidomain coupling manager, which will couple the mass and the momentum problems +We can obtain the type from either the `MomentumTypeTag` or the `MassTypeTag` because they are mutually coupled with the same manager ```cpp using CouplingManager = GetPropType; diff --git a/test/discretization/cellcentered/tpfa/test_tpfafvgeometry_nonconforming.cc b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry_nonconforming.cc index 0b5b2c380db7a61ccf4640010f64dc53c61167a4..de1c1d61a612a39d386caa9e221264af2bfdf14c 100644 --- a/test/discretization/cellcentered/tpfa/test_tpfafvgeometry_nonconforming.cc +++ b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry_nonconforming.cc @@ -203,8 +203,8 @@ int main (int argc, char *argv[]) DUNE_THROW(Dune::InvalidStateException, "Scvf has more than one neighbor"); //! center must always be between the corners - const auto d1 = scvf.corner(0) - scvf.center(); - const auto d2 = scvf.corner(1) - scvf.center(); + const auto d1 = fvGeometry.geometry(scvf).corner(0) - scvf.center(); + const auto d2 = fvGeometry.geometry(scvf).corner(1) - scvf.center(); if ( d1 * d2 >= 0 ) DUNE_THROW(Dune::InvalidStateException, "Center is not between the two corners"); diff --git a/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc b/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc index c761e86fa5e974cda207ac64a1a1f4925d529fbf..fa54ce4cadc6904041554e5945b439dd1daec5ef 100644 --- a/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc +++ b/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc @@ -48,35 +48,35 @@ void runTest(const GG& gg, const double refVolume, const double refSurface) fvGeometry.bind(element); for (const auto& scv : scvs(fvGeometry)) - volume += GG::Extrusion::volume(scv); + volume += GG::Extrusion::volume(fvGeometry, scv); for (const auto& scvf : scvfs(fvGeometry)) if (scvf.boundary()) - surface += GG::Extrusion::area(scvf); + surface += GG::Extrusion::area(fvGeometry, scvf); // compare volume and integrated volume of one scv const auto& scv = *(scvs(fvGeometry).begin()); - const auto scvGeometry = scv.geometry(); + const auto scvGeometry = fvGeometry.geometry(scv); double volScv = 0.0; const auto ruleScv = Dune::QuadratureRules::rule(scvGeometry.type(), 3); for (const auto& qp : ruleScv) volScv += qp.weight()*GG::Extrusion::integrationElement(scvGeometry, qp.position()); - if (!Dune::FloatCmp::eq(volScv, GG::Extrusion::volume(scv))) - DUNE_THROW(Dune::Exception, "Integration not correct! Integrated: " << volScv << ", direct: " << GG::Extrusion::volume(scv)); + if (!Dune::FloatCmp::eq(volScv, GG::Extrusion::volume(fvGeometry, scv))) + DUNE_THROW(Dune::Exception, "Integration not correct! Integrated: " << volScv << ", direct: " << GG::Extrusion::volume(fvGeometry, scv)); // compare area and integration area of one scvf const auto& scvf = *(scvfs(fvGeometry).begin()); - const auto scvfGeometry = scvf.geometry(); + const auto scvfGeometry = fvGeometry.geometry(scvf); double volScvf = 0.0; const auto ruleScvf = Dune::QuadratureRules::rule(scvfGeometry.type(), 3); for (const auto& qp : ruleScvf) volScvf += qp.weight()*GG::Extrusion::integrationElement(scvfGeometry, qp.position()); - if (!Dune::FloatCmp::eq(volScvf, GG::Extrusion::area(scvf))) - DUNE_THROW(Dune::Exception, "Integration not correct! Integrated: " << volScvf << ", direct: " << GG::Extrusion::area(scvf)); + if (!Dune::FloatCmp::eq(volScvf, GG::Extrusion::area(fvGeometry, scvf))) + DUNE_THROW(Dune::Exception, "Integration not correct! Integrated: " << volScvf << ", direct: " << GG::Extrusion::area(fvGeometry, scvf)); } // compare total volume/surface area to reference diff --git a/test/freeflow/navierstokes/errors.hh b/test/freeflow/navierstokes/errors.hh index 89e872a83a18b759e29ba2e5b90d953b4485c1ef..f8a0295838a9a024ca567802b23a5771016e3ea3 100644 --- a/test/freeflow/navierstokes/errors.hh +++ b/test/freeflow/navierstokes/errors.hh @@ -102,7 +102,7 @@ private: for (const auto& scv : scvs(fvGeometry)) { - totalVolume += Extrusion::volume(scv); + totalVolume += Extrusion::volume(fvGeometry, scv); // compute the pressure errors const auto analyticalSolutionCellCenter @@ -115,17 +115,17 @@ private: maxError[Indices::pressureIdx] = std::max(maxError[Indices::pressureIdx], pError); maxReference[Indices::pressureIdx] = std::max(maxReference[Indices::pressureIdx], pReference); - sumError[Indices::pressureIdx] += pError * pError * Extrusion::volume(scv); - sumReference[Indices::pressureIdx] += pReference * pReference * Extrusion::volume(scv); + sumError[Indices::pressureIdx] += pError * pError * Extrusion::volume(fvGeometry, scv); + sumReference[Indices::pressureIdx] += pReference * pReference * Extrusion::volume(fvGeometry, scv); for (const auto& scvf : scvfs(fvGeometry)) { // compute the velocity errors const auto velocityIndex = Indices::velocity(scvf.directionIndex()); - const Scalar staggeredHalfVolume = Extrusion::volume( + const Scalar staggeredHalfVolume = Extrusion::volume(fvGeometry, typename GridGeometry::Traits::FaceSubControlVolume( - 0.5*(scv.center() + scvf.center()), 0.5*Extrusion::volume(scv) + 0.5*(scv.center() + scvf.center()), 0.5*Extrusion::volume(fvGeometry, scv) ) ); @@ -376,7 +376,7 @@ private: { using GridGeometry = std::decay_t().gridGeometry())>; using Extrusion = Extrusion_t; - totalVolume_ += Extrusion::volume(scv); + totalVolume_ += Extrusion::volume(fvGeometry, scv); // velocity errors if constexpr (isMomentumProblem) @@ -396,8 +396,8 @@ private: maxError[velIdx] = std::max(maxError[velIdx], vError); maxReference[velIdx] = std::max(maxReference[velIdx], vReference); - sumError[velIdx] += vError * vError * Extrusion::volume(scv); - sumReference[velIdx] += vReference * vReference * Extrusion::volume(scv); + sumError[velIdx] += vError * vError * Extrusion::volume(fvGeometry, scv); + sumReference[velIdx] += vReference * vReference * Extrusion::volume(fvGeometry, scv); } else if (GridGeometry::discMethod == DiscretizationMethods::fcdiamond) { @@ -413,8 +413,8 @@ private: maxError[dirIdx] = std::max(maxError[dirIdx], vError); maxReference[dirIdx] = std::max(maxReference[dirIdx], vReference); - sumError[dirIdx] += vError * vError * Extrusion::volume(scv); - sumReference[dirIdx] += vReference * vReference * Extrusion::volume(scv); + sumError[dirIdx] += vError * vError * Extrusion::volume(fvGeometry, scv); + sumReference[dirIdx] += vReference * vReference * Extrusion::volume(fvGeometry, scv); } } } @@ -433,8 +433,8 @@ private: maxError[0] = std::max(maxError[0], pError); maxReference[0] = std::max(maxReference[0], pReference); - sumError[0] += pError * pError * Extrusion::volume(scv); - sumReference[0] += pReference * pReference * Extrusion::volume(scv); + sumError[0] += pError * pError * Extrusion::volume(fvGeometry, scv); + sumReference[0] += pReference * pReference * Extrusion::volume(fvGeometry, scv); } } } diff --git a/test/multidomain/embedded/1d3d/1p_1p/problem_bloodflow.hh b/test/multidomain/embedded/1d3d/1p_1p/problem_bloodflow.hh index 6632b0b12402decd127687b6b05daac338405aaf..d82ff87310a93155be54ef769b5c5d9eb82e5868 100644 --- a/test/multidomain/embedded/1d3d/1p_1p/problem_bloodflow.hh +++ b/test/multidomain/embedded/1d3d/1p_1p/problem_bloodflow.hh @@ -225,8 +225,8 @@ public: auto pointSources = this->scvPointSources(element, fvGeometry, elemVolVars, scv); pointSources *= scv.volume()*elemVolVars[scv].extrusionFactor(); source[scv.dofIndex()] += pointSources; - auto a = scv.corner(0)[2]; - auto b = scv.corner(1)[2]; + auto a = fvGeometry.geometry(scv).corner(0)[2]; + auto b = fvGeometry.geometry(scv).corner(1)[2]; if (a > b) std::swap(a, b); sourceExact[scv.dofIndex()] += a*(1.0+0.5*a) - b*(1.0+0.5*b); volume += scv.volume(); diff --git a/test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh b/test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh index 0402757d710a1710dc56deb6af12adcc919294bb..c8833ff21434c3bd2c4bbc0263b4c547b3982c57 100644 --- a/test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh +++ b/test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh @@ -157,7 +157,7 @@ public: { NumEqVector flux(0.0); // integrate over the scvf to compute the flux - const auto geometry = scvf.geometry(); + const auto geometry = fvGeometry.geometry(scvf); Scalar derivative = 0.0; const auto& quad = Dune::QuadratureRules::rule(geometry.type(), 4); for(auto&& qp : quad) diff --git a/test/porousmediumflow/richards/annulus/problem.hh b/test/porousmediumflow/richards/annulus/problem.hh index 1996c5f8f419aec58e7b2f7dbbb9f13c1ad11166..2104e254109255dfbe58ed786731ae47c252cca5 100644 --- a/test/porousmediumflow/richards/annulus/problem.hh +++ b/test/porousmediumflow/richards/annulus/problem.hh @@ -194,8 +194,8 @@ public: for (const auto& scv : scvs(fvGeometry)) { const auto& volVars = elemVolVars[scv]; - totalWaterVolume += Extrusion::volume(scv)*volVars.porosity()*volVars.saturation(0); - refWaterVolume += Extrusion::volume(scv)*volVars.porosity()*initialSaturation; + totalWaterVolume += Extrusion::volume(fvGeometry, scv)*volVars.porosity()*volVars.saturation(0); + refWaterVolume += Extrusion::volume(fvGeometry, scv)*volVars.porosity()*initialSaturation; } } @@ -223,8 +223,8 @@ public: { const auto& volVars = elemVolVars[scv]; const auto& oldVolVars = oldElemVolVars[scv]; - storageDerivative[scv.dofIndex()] += Extrusion::volume(scv)*(volVars.saturation(0) - oldVolVars.saturation(0))/dt; - volumes[scv.dofIndex()] += Extrusion::volume(scv); + storageDerivative[scv.dofIndex()] += Extrusion::volume(fvGeometry, scv)*(volVars.saturation(0) - oldVolVars.saturation(0))/dt; + volumes[scv.dofIndex()] += Extrusion::volume(fvGeometry, scv); } }