From 16550a66e5b3a751e274905853baea69ef53aa89 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 14 Sep 2022 10:37:52 +0000 Subject: [PATCH 01/20] [geometry] Add volume function with transformed integrationElement This can be used e.g. if there is an extrusion transformation involved --- dumux/geometry/volume.hh | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/dumux/geometry/volume.hh b/dumux/geometry/volume.hh index 9650d2ca8b..4883cf8010 100644 --- a/dumux/geometry/volume.hh +++ b/dumux/geometry/volume.hh @@ -177,6 +177,21 @@ auto volume(const Geometry& geo, unsigned int integrationOrder = 4) 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) +{ + double 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 -- GitLab From b61c06329cf8fc8e847c5189db10636ec551bf28 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 14 Sep 2022 10:41:28 +0000 Subject: [PATCH 02/20] [geometry] Deduce coordinate type in volume --- dumux/geometry/volume.hh | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/dumux/geometry/volume.hh b/dumux/geometry/volume.hh index 4883cf8010..ea2888714f 100644 --- a/dumux/geometry/volume.hh +++ b/dumux/geometry/volume.hh @@ -170,8 +170,9 @@ 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; @@ -185,8 +186,9 @@ auto volume(const Geometry& geo, unsigned int integrationOrder = 4) template auto volume(const Geometry& geo, Transformation transformation, 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 += transformation.integrationElement(geo, qp.position())*qp.weight(); return volume; -- GitLab From 99d6445d5778435d981214efb18be6efc9f6dce2 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Tue, 13 Sep 2022 09:19:30 +0200 Subject: [PATCH 03/20] Use new area/volume with fvGeometry argument --- dumux/assembly/cclocalresidual.hh | 2 +- dumux/assembly/cvfelocalresidual.hh | 4 ++-- dumux/assembly/fclocalresidual.hh | 2 +- dumux/assembly/fvlocalresidual.hh | 6 ++--- dumux/assembly/staggeredlocalresidual.hh | 4 ++-- dumux/common/fvproblem.hh | 2 +- .../cellcentered/mpfa/localassemblerbase.hh | 6 ++--- dumux/flux/box/darcyslaw.hh | 4 ++-- dumux/flux/box/dispersionflux.hh | 2 +- dumux/flux/box/effectivestresslaw.hh | 2 +- dumux/flux/box/fickslaw.hh | 2 +- dumux/flux/box/forchheimerslaw.hh | 2 +- dumux/flux/box/fourierslaw.hh | 2 +- dumux/flux/box/fourierslawnonequilibrium.hh | 2 +- dumux/flux/box/hookeslaw.hh | 2 +- dumux/flux/box/maxwellstefanslaw.hh | 2 +- dumux/flux/cctpfa/darcyslaw.hh | 16 +++++++------- dumux/flux/cctpfa/fickslaw.hh | 2 +- dumux/flux/cctpfa/forchheimerslaw.hh | 4 ++-- dumux/flux/cctpfa/fourierslaw.hh | 4 ++-- .../flux/cctpfa/fourierslawnonequilibrium.hh | 4 ++-- dumux/flux/cctpfa/maxwellstefanslaw.hh | 2 +- dumux/flux/cvfe/darcyslaw.hh | 4 ++-- dumux/flux/staggered/freeflow/fickslaw.hh | 2 +- dumux/flux/staggered/freeflow/fourierslaw.hh | 2 +- .../staggered/freeflow/maxwellstefanslaw.hh | 2 +- .../navierstokes/momentum/diamond/flux.hh | 6 ++--- .../navierstokes/momentum/fluxvariables.hh | 12 +++++----- .../navierstokes/momentum/localresidual.hh | 6 ++--- .../navierstokes/scalarfluxvariables.hh | 4 ++-- .../navierstokes/staggered/fluxoversurface.hh | 2 +- .../navierstokes/staggered/localresidual.hh | 6 ++--- .../rans/oneeq/staggered/fluxvariables.hh | 2 +- .../twoeq/kepsilon/staggered/fluxvariables.hh | 6 ++--- .../twoeq/komega/staggered/fluxvariables.hh | 6 ++--- .../lowrekepsilon/staggered/fluxvariables.hh | 6 ++--- .../rans/twoeq/sst/staggered/fluxvariables.hh | 6 ++--- .../boundary/darcydarcy/couplingmanager.hh | 6 ++--- .../embedded/couplingmanagerbase.hh | 2 +- dumux/multidomain/facet/box/darcyslaw.hh | 4 ++-- dumux/multidomain/facet/box/fickslaw.hh | 4 ++-- dumux/multidomain/facet/box/fourierslaw.hh | 4 ++-- dumux/multidomain/facet/box/localresidual.hh | 2 +- .../facet/cellcentered/mpfa/localassembler.hh | 8 +++---- .../facet/cellcentered/tpfa/darcyslaw.hh | 12 +++++----- .../facet/cellcentered/tpfa/fickslaw.hh | 10 ++++----- .../facet/cellcentered/tpfa/fourierslaw.hh | 10 ++++----- .../multidomain/subdomainboxlocalassembler.hh | 2 +- .../multidomain/subdomaincclocalassembler.hh | 2 +- dumux/porousmediumflow/2p/griddatatransfer.hh | 12 +++++----- .../2p/incompressiblelocalresidual.hh | 2 +- .../fluxvariablescachefiller.hh | 2 +- .../richards/localresidual.hh | 2 +- .../porousmediumflow/tracer/localresidual.hh | 2 +- dumux/porousmediumflow/velocity.hh | 4 ++-- .../test_rotationsymmetric_gridgeometry.cc | 12 +++++----- test/freeflow/navierstokes/errors.hh | 22 +++++++++---------- .../richards/annulus/problem.hh | 8 +++---- 58 files changed, 141 insertions(+), 141 deletions(-) diff --git a/dumux/assembly/cclocalresidual.hh b/dumux/assembly/cclocalresidual.hh index e69b73f760..e02236c6b9 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 00122079f3..a360aa32e5 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 5d4ff9051f..15a2cca28e 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 55872ff55e..b5d0fea004 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 927f01d144..f13771715b 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; diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index bf040be778..a814d2a49c 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/localassemblerbase.hh b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh index af558765ae..1e8d13f80a 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/flux/box/darcyslaw.hh b/dumux/flux/box/darcyslaw.hh index e01000f1b3..ceb511b949 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 b1a6c6d1ca..0a3b7dfe81 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 ebd8d4ac59..de8271e9e4 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 46c5badcf8..394b15ca27 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; diff --git a/dumux/flux/box/forchheimerslaw.hh b/dumux/flux/box/forchheimerslaw.hh index f7c1288642..ae6e3decb1 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 e91d6503a3..f29d3b9460 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 1c5024088b..f0118f6b82 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 f750acd26f..d29dfe167d 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 4e21b7eb17..cdea991878 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 fa9b7eb164..88a96791c2 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()) @@ -245,7 +245,7 @@ class CCTpfaDarcysLaw // 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 @@ -264,7 +264,7 @@ class CCTpfaDarcysLaw 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) @@ -487,7 +487,7 @@ public: // 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 @@ -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/fickslaw.hh b/dumux/flux/cctpfa/fickslaw.hh index ce5847b4b8..970dcd4f6c 100644 --- a/dumux/flux/cctpfa/fickslaw.hh +++ b/dumux/flux/cctpfa/fickslaw.hh @@ -246,7 +246,7 @@ public: // 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 diff --git a/dumux/flux/cctpfa/forchheimerslaw.hh b/dumux/flux/cctpfa/forchheimerslaw.hh index 7e1065254d..5e05afe757 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 @@ -209,7 +209,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/fourierslawnonequilibrium.hh b/dumux/flux/cctpfa/fourierslawnonequilibrium.hh index a2ddb1e83d..ddd7fbde5b 100644 --- a/dumux/flux/cctpfa/fourierslawnonequilibrium.hh +++ b/dumux/flux/cctpfa/fourierslawnonequilibrium.hh @@ -129,7 +129,7 @@ public: // 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(); @@ -148,7 +148,7 @@ public: 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 c45578cb0c..9e08cb09eb 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 42ec2a1047..ed05f5c297 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 96359c313d..ae4a8b1885 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 1221ec07c1..76d2e2230c 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 d6408f7893..66960d446e 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/navierstokes/momentum/diamond/flux.hh b/dumux/freeflow/navierstokes/momentum/diamond/flux.hh index e1b47b96f7..03c3421c9f 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 4d350f9e3e..6cc1f97c56 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 43924b431f..68cb238c38 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 3a07b5e9bb..f440ca34ba 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 219bab9561..3873fc1564 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/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh index 55d905da70..7c363236cc 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()); @@ -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/rans/oneeq/staggered/fluxvariables.hh b/dumux/freeflow/rans/oneeq/staggered/fluxvariables.hh index 8b18e0d18a..246292b05f 100644 --- a/dumux/freeflow/rans/oneeq/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/oneeq/staggered/fluxvariables.hh @@ -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 a7232f9e05..1988732e92 100644 --- a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh @@ -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 e8b86e24b2..ced76ff43e 100644 --- a/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh @@ -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 b6aa6f9086..3e49ebb929 100644 --- a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh @@ -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 d6f62f6146..70b17f9098 100644 --- a/dumux/freeflow/rans/twoeq/sst/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/twoeq/sst/staggered/fluxvariables.hh @@ -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/multidomain/boundary/darcydarcy/couplingmanager.hh b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh index 8610540439..09f845989b 100644 --- a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh +++ b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh @@ -182,7 +182,7 @@ public: const auto tj = computeTpfaTransmissibility(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/embedded/couplingmanagerbase.hh b/dumux/multidomain/embedded/couplingmanagerbase.hh index ff47938755..b691cbac07 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 45b460a338..d404520cf6 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 075aa5f175..118107a84b 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 a6363228c3..5bc2432b01 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 7d1d2568e9..cfba0bbcfc 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/cellcentered/mpfa/localassembler.hh b/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh index d28b08cb6b..7c063771c8 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); } } @@ -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()); diff --git a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh index 05c4ec5d7a..42c9691ef7 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 fb79b537fc..5cbb38c64b 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/porousmediumflow/2p/griddatatransfer.hh b/dumux/porousmediumflow/2p/griddatatransfer.hh index 38ec15039d..5dc596abe1 100644 --- a/dumux/porousmediumflow/2p/griddatatransfer.hh +++ b/dumux/porousmediumflow/2p/griddatatransfer.hh @@ -163,7 +163,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); } @@ -272,7 +272,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 +324,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)/Extrusion::volume(fatherElement.geometry())*massFather )/massCoeffSon; } } else @@ -354,7 +354,7 @@ public: 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 659f70db00..922fed8670 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/fluxvariablescachefiller.hh b/dumux/porousmediumflow/fluxvariablescachefiller.hh index 1ec42a8d13..ff2d633142 100644 --- a/dumux/porousmediumflow/fluxvariablescachefiller.hh +++ b/dumux/porousmediumflow/fluxvariablescachefiller.hh @@ -637,7 +637,7 @@ 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( + const auto eps = Extrusion::area(fvGeometry(), scvf)*computeTpfaTransmissibility( scvf, scv, zeroD, vv.extrusionFactor() ); diff --git a/dumux/porousmediumflow/richards/localresidual.hh b/dumux/porousmediumflow/richards/localresidual.hh index 218c06a5be..a3fc5ed588 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 754d5b568c..770674990e 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 8f0af5f702..6f89955df4 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/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc b/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc index c761e86fa5..14771352d4 100644 --- a/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc +++ b/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc @@ -48,11 +48,11 @@ 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()); @@ -63,8 +63,8 @@ void runTest(const GG& gg, const double refVolume, const double refSurface) 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()); @@ -75,8 +75,8 @@ void runTest(const GG& gg, const double refVolume, const double refSurface) 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 89e872a83a..ddc5925aa0 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,8 +115,8 @@ 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)) { @@ -125,7 +125,7 @@ private: const Scalar staggeredHalfVolume = Extrusion::volume( 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/porousmediumflow/richards/annulus/problem.hh b/test/porousmediumflow/richards/annulus/problem.hh index 1996c5f8f4..2104e25410 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); } } -- GitLab From cf370eeebce5839508dd5279b54cc9526f662fd0 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Tue, 13 Sep 2022 09:21:24 +0200 Subject: [PATCH 04/20] [disc] Use second template param in Extrusion::area/volume --- dumux/discretization/extrusion.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dumux/discretization/extrusion.hh b/dumux/discretization/extrusion.hh index 19fabb9e85..8e592f3e01 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& fvGeometry, 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& fvGeometry, const SCV& scv) { return scv.volume(); } template -- GitLab From 8941ac814ae1bbedb1277e1575749a82bb762465 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Tue, 13 Sep 2022 09:22:06 +0200 Subject: [PATCH 05/20] [disc] Add fvGeometry.geometry(scv/scvf) interface --- .../cellcentered/tpfa/fvelementgeometry.hh | 24 +++++++++++++++++++ .../staggered/fvelementgeometry.hh | 6 +++++ .../staggered/fvelementgeometry.hh | 8 ++++++- 3 files changed, 37 insertions(+), 1 deletion(-) diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh index 6129f321aa..e36175f2f4 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/facecentered/staggered/fvelementgeometry.hh b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh index 9f08c7e6b1..2f62de3574 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 6581a68b9c..0f37a11550 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) -- GitLab From 7f34a0d1dd891d81870d427239b6d003d39cd3ee Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Tue, 13 Sep 2022 09:23:00 +0200 Subject: [PATCH 06/20] [porenetwork][2p] Use cvfe headers for fluxvariablescaches --- dumux/porenetwork/2p/elementfluxvariablescache.hh | 6 +++--- dumux/porenetwork/2p/gridfluxvariablescache.hh | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/dumux/porenetwork/2p/elementfluxvariablescache.hh b/dumux/porenetwork/2p/elementfluxvariablescache.hh index e111fc1eb9..92bd6b8bbd 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 6b05046147..1dcb53ac5c 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" -- GitLab From 1ccb10d812dc4ac0690629f8cc9990ff36c60f88 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Tue, 13 Sep 2022 09:23:30 +0200 Subject: [PATCH 07/20] Use new Extrusion/fvGeometry interfaces --- dumux/assembly/staggeredlocalresidual.hh | 4 ++-- dumux/discretization/walldistance.hh | 2 +- dumux/flux/box/fickslaw.hh | 2 +- dumux/flux/cctpfa/fickslaw.hh | 2 +- .../navierstokes/staggered/fluxvariables.hh | 18 +++++++++--------- .../freeflowporenetwork/couplingmapper.hh | 2 +- .../ffmomentumpm/couplingmapper.hh | 4 ++-- .../tpfa/test_tpfafvgeometry_nonconforming.cc | 4 ++-- .../test_rotationsymmetric_gridgeometry.cc | 4 ++-- test/freeflow/navierstokes/errors.hh | 2 +- .../embedded/1d3d/1p_1p/problem_bloodflow.hh | 4 ++-- .../embedded/1d3d/1p_1p/problem_tissue.hh | 2 +- 12 files changed, 25 insertions(+), 25 deletions(-) diff --git a/dumux/assembly/staggeredlocalresidual.hh b/dumux/assembly/staggeredlocalresidual.hh index f13771715b..15826b6bbf 100644 --- a/dumux/assembly/staggeredlocalresidual.hh +++ b/dumux/assembly/staggeredlocalresidual.hh @@ -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/discretization/walldistance.hh b/dumux/discretization/walldistance.hh index 045b8e2ce6..ae485bee92 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/fickslaw.hh b/dumux/flux/box/fickslaw.hh index 394b15ca27..f94549f3ee 100644 --- a/dumux/flux/box/fickslaw.hh +++ b/dumux/flux/box/fickslaw.hh @@ -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/cctpfa/fickslaw.hh b/dumux/flux/cctpfa/fickslaw.hh index 970dcd4f6c..1a65750bd9 100644 --- a/dumux/flux/cctpfa/fickslaw.hh +++ b/dumux/flux/cctpfa/fickslaw.hh @@ -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/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh index 2ab1a75a9d..1af5932268 100644 --- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh @@ -130,7 +130,7 @@ public: const Scalar flux = (upwindWeight * upwindTerm(upstreamVolVars) + (1.0 - upwindWeight) * upwindTerm(downstreamVolVars)) - * velocity * Extrusion::area(scvf) * scvf.directionSign(); + * velocity * Extrusion::area(localView(problem.gridGeometry()), scvf) * scvf.directionSign(); return flux * extrusionFactor_(elemVolVars, scvf); } @@ -255,7 +255,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 +334,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 +350,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 +374,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; } @@ -459,7 +459,7 @@ public: inOrOutflow += boundaryPressure; // Account for the orientation of the face at the boundary, - return inOrOutflow * scvf.directionSign() * Extrusion::area(scvf) * insideVolVars.extrusionFactor(); + return inOrOutflow * scvf.directionSign() * Extrusion::area(localView(problem.gridGeometry()), scvf) * insideVolVars.extrusionFactor(); } private: @@ -533,7 +533,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 +611,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 +688,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/multidomain/boundary/freeflowporenetwork/couplingmapper.hh b/dumux/multidomain/boundary/freeflowporenetwork/couplingmapper.hh index c8d5e17e23..6b10c18b71 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/ffmomentumpm/couplingmapper.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper.hh index 5bca370ef4..1e64decf02 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/test/discretization/cellcentered/tpfa/test_tpfafvgeometry_nonconforming.cc b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry_nonconforming.cc index 0b5b2c380d..de1c1d61a6 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 14771352d4..fa54ce4cad 100644 --- a/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc +++ b/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc @@ -56,7 +56,7 @@ void runTest(const GG& gg, const double refVolume, const double refSurface) // 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); @@ -68,7 +68,7 @@ void runTest(const GG& gg, const double refVolume, const double refSurface) // 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); diff --git a/test/freeflow/navierstokes/errors.hh b/test/freeflow/navierstokes/errors.hh index ddc5925aa0..f8a0295838 100644 --- a/test/freeflow/navierstokes/errors.hh +++ b/test/freeflow/navierstokes/errors.hh @@ -123,7 +123,7 @@ private: // 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(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 6632b0b124..d82ff87310 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 0402757d71..c8833ff214 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) -- GitLab From a9468879d3739156f8f66064bbec527fa7b864d1 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Tue, 13 Sep 2022 09:24:01 +0200 Subject: [PATCH 08/20] [porousmediumflow] Use geometry.volume instead of Extrusion --- dumux/porousmediumflow/2p/griddatatransfer.hh | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/dumux/porousmediumflow/2p/griddatatransfer.hh b/dumux/porousmediumflow/2p/griddatatransfer.hh index 5dc596abe1..43d4d33547 100644 --- a/dumux/porousmediumflow/2p/griddatatransfer.hh +++ b/dumux/porousmediumflow/2p/griddatatransfer.hh @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -240,7 +241,7 @@ public: if (!isBox) elemSol[0] /= adaptedValues.count; - const auto elementVolume = Extrusion::volume(element.geometry()); + const auto elementVolume = volume(element.geometry()); for (const auto& scv : scvs(fvGeometry)) { VolumeVariables volVars; @@ -328,7 +329,7 @@ public: else if (formulation == p1s0) massCoeffSon = Extrusion::volume(fvGeometry, scv) * volVars.density(phase0Idx) * volVars.porosity(); sol_[scv.dofIndex()][saturationIdx] = - ( Extrusion::volume(fvGeometry, scv)/Extrusion::volume(fatherElement.geometry())*massFather )/massCoeffSon; + ( Extrusion::volume(fvGeometry, scv)/volume(fatherElement.geometry())*massFather )/massCoeffSon; } } else @@ -347,7 +348,7 @@ 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); for (const auto& scv : scvs(fvGeometry)) { VolumeVariables volVars; -- GitLab From 766db84d3f665d7247c2e634373a1680669d0b3c Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Tue, 13 Sep 2022 09:24:28 +0200 Subject: [PATCH 09/20] temp remaining deprecation warnings --- .../cellcentered/mpfa/computetransmissibility.hh | 4 ++-- dumux/multidomain/facet/box/subcontrolvolumeface.hh | 4 ++-- dumux/porousmediumflow/boxdfm/subcontrolvolume.hh | 4 ++-- dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh b/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh index 0159500ac7..c1bb46d0e3 100644 --- a/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh +++ b/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh @@ -58,7 +58,7 @@ computeMpfaTransmissibility(const IVSubControlVolume& scv, wijk[dir] = vtmv(scvf.unitOuterNormal(), t, scv.nu(dir)); using Extrusion = Extrusion_t; - wijk *= Extrusion::area(scvf)*extrusionFactor; + wijk *= Extrusion::area(scvf)*extrusionFactor;// TODO: No fvGeometry available wijk /= scv.detX(); return wijk; @@ -88,7 +88,7 @@ computeMpfaTransmissibility(const IVSubControlVolume& scv, wijk[dir] = vtmv(scvf.unitOuterNormal(), t, scv.nu(dir)); using Extrusion = Extrusion_t; - wijk *= Extrusion::area(scvf)*extrusionFactor; + wijk *= Extrusion::area(scvf)*extrusionFactor;// TODO: No fvGeometry available wijk /= scv.detX(); return wijk; diff --git a/dumux/multidomain/facet/box/subcontrolvolumeface.hh b/dumux/multidomain/facet/box/subcontrolvolumeface.hh index 207de2a27b..aa023ff090 100644 --- a/dumux/multidomain/facet/box/subcontrolvolumeface.hh +++ b/dumux/multidomain/facet/box/subcontrolvolumeface.hh @@ -78,7 +78,7 @@ public: : corners_(geometryHelper.getScvfCorners(scvfIndex)) , center_(0.0) , unitOuterNormal_(geometryHelper.normal(corners_, scvIndices)) - , area_(geometryHelper.scvfArea(corners_)) + , area_(geometryHelper.scvfArea(corners_))//TODO , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) , facetIndex_(/*undefined*/) @@ -105,7 +105,7 @@ public: : corners_(geometryHelper.getBoundaryScvfCorners(intersection.indexInInside(), indexInIntersection)) , center_(0.0) , unitOuterNormal_(intersection.centerUnitOuterNormal()) - , area_(geometryHelper.scvfArea(corners_)) + , area_(geometryHelper.scvfArea(corners_))//TODO , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) , facetIndex_(intersection.indexInInside()) diff --git a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh index cb880d192c..b06e1c676d 100644 --- a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh +++ b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh @@ -101,7 +101,7 @@ public: : isFractureScv_(false) , corners_(geometryHelper.getScvCorners(scvIdx)) , center_(0.0) - , volume_(geometryHelper.scvVolume(corners_)) + , volume_(geometryHelper.scvVolume(corners_))//TODO , elementIndex_(elementIndex) , vIdxLocal_(scvIdx) , elemLocalScvIdx_(scvIdx) @@ -151,7 +151,7 @@ public: corners_[i] = corners[i]; // compute volume and scv center - volume_ = geometryHelper.scvfArea(corners); + volume_ = geometryHelper.scvfArea(corners);//TODO 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 6c7f0df124..a404c7826e 100644 --- a/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh +++ b/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh @@ -103,7 +103,7 @@ public: : corners_(geometryHelper.getScvfCorners(scvfIndex)) , center_(0.0) , unitOuterNormal_(geometryHelper.normal(corners_, scvIndices)) - , area_(geometryHelper.scvfArea(corners_)) + , area_(geometryHelper.scvfArea(corners_))//TODO , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) , boundary_(false) @@ -127,7 +127,7 @@ public: : corners_(geometryHelper.getBoundaryScvfCorners(intersection.indexInInside(), indexInIntersection)) , center_(0.0) , unitOuterNormal_(intersection.centerUnitOuterNormal()) - , area_(geometryHelper.scvfArea(corners_)) + , area_(geometryHelper.scvfArea(corners_))//TODO , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) , boundary_(true) -- GitLab From da67443d76731ecfb3c9304c379be27eeec91ff0 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Tue, 13 Sep 2022 13:48:29 +0200 Subject: [PATCH 10/20] [examples] fixup adjust documentation --- examples/liddrivencavity/doc/problem.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/liddrivencavity/doc/problem.md b/examples/liddrivencavity/doc/problem.md index e70cd12f88..e2d607cd1e 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; -- GitLab From 0789c0aaf5860e03e447335837f62986711ea4a2 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Tue, 13 Sep 2022 18:17:30 +0200 Subject: [PATCH 11/20] [box] Use convexPolytopeVolume in scv(f) constructors --- dumux/multidomain/facet/box/subcontrolvolumeface.hh | 11 +++++++++-- dumux/porousmediumflow/boxdfm/subcontrolvolume.hh | 11 +++++++++-- dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh | 11 +++++++++-- 3 files changed, 27 insertions(+), 6 deletions(-) diff --git a/dumux/multidomain/facet/box/subcontrolvolumeface.hh b/dumux/multidomain/facet/box/subcontrolvolumeface.hh index aa023ff090..0c024c883c 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_))//TODO + , 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_))//TODO + , 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/porousmediumflow/boxdfm/subcontrolvolume.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh index b06e1c676d..637932dfe7 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_))//TODO + , 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);//TODO + 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 a404c7826e..8cc86bf696 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_))//TODO + , 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_))//TODO + , area_(Dumux::convexPolytopeVolume( + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; }) + ) , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) , boundary_(true) -- GitLab From 034cddfce815d7582d39ed464f97f84cedb17056 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Wed, 14 Sep 2022 09:29:23 +0200 Subject: [PATCH 12/20] [freeflow] add fvGeometry parameter to staggered fluxvars interfaces --- .../compositional/staggered/fluxvariables.hh | 2 +- .../navierstokes/staggered/fluxvariables.hh | 89 ++++++++++++++++++- .../navierstokes/staggered/localresidual.hh | 2 +- dumux/freeflow/nonisothermal/localresidual.hh | 1 + .../rans/oneeq/staggered/fluxvariables.hh | 2 +- .../twoeq/kepsilon/staggered/fluxvariables.hh | 4 +- .../twoeq/komega/staggered/fluxvariables.hh | 4 +- .../lowrekepsilon/staggered/fluxvariables.hh | 4 +- .../rans/twoeq/sst/staggered/fluxvariables.hh | 4 +- 9 files changed, 98 insertions(+), 14 deletions(-) diff --git a/dumux/freeflow/compositional/staggered/fluxvariables.hh b/dumux/freeflow/compositional/staggered/fluxvariables.hh index b6c3db8f51..cf9bf0659f 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/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh index 1af5932268..0945f55fcd 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, @@ -130,7 +131,40 @@ public: const Scalar flux = (upwindWeight * upwindTerm(upstreamVolVars) + (1.0 - upwindWeight) * upwindTerm(downstreamVolVars)) - * velocity * Extrusion::area(localView(problem.gridGeometry()), scvf) * scvf.directionSign(); + * velocity * Extrusion::area(scvf) * scvf.directionSign(); + + 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 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); } @@ -163,7 +197,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; } @@ -415,6 +449,54 @@ public: return lateralFlux; } + /*! + * \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 + */ + [[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, + const ElementVolumeVariables& elemVolVars, + const ElementFaceVariables& elemFaceVars) const + { + FacePrimaryVariables inOrOutflow(0.0); + 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(scvf) * insideVolVars.extrusionFactor(); + } + /*! * \brief Returns the momentum flux over an inflow or outflow boundary face. * @@ -434,6 +516,7 @@ public: FacePrimaryVariables inflowOutflowBoundaryFlux(const Problem& problem, const Element& element, const SubControlVolumeFace& scvf, + const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const ElementFaceVariables& elemFaceVars) const { @@ -459,7 +542,7 @@ public: inOrOutflow += boundaryPressure; // Account for the orientation of the face at the boundary, - return inOrOutflow * scvf.directionSign() * Extrusion::area(localView(problem.gridGeometry()), scvf) * insideVolVars.extrusionFactor(); + return inOrOutflow * scvf.directionSign() * Extrusion::area(fvGeometry, scvf) * insideVolVars.extrusionFactor(); } private: diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh index 7c363236cc..04f85e3008 100644 --- a/dumux/freeflow/navierstokes/staggered/localresidual.hh +++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh @@ -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, element, scvf, fvGeometry, elemVolVars, elemFaceVars); } } return result; diff --git a/dumux/freeflow/nonisothermal/localresidual.hh b/dumux/freeflow/nonisothermal/localresidual.hh index 2a5bc519a9..d972ae0084 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 246292b05f..0c4289a521 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()); diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh index 1988732e92..2a297d3a4b 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()); diff --git a/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh index ced76ff43e..e347aa55d7 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()); diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh index 3e49ebb929..874dd54aef 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()); diff --git a/dumux/freeflow/rans/twoeq/sst/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/sst/staggered/fluxvariables.hh index 70b17f9098..8384e2b485 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()); -- GitLab From 9d3e3d188c5021ebf0ea45aba7e77aaf43ab6174 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Wed, 14 Sep 2022 14:10:10 +0200 Subject: [PATCH 13/20] [disc] Add fvGeometry parameter to transmissibilities --- .../mpfa/computetransmissibility.hh | 68 ++++++++++++++++++- .../mpfa/omethod/localassembler.hh | 6 +- .../tpfa/computetransmissibility.hh | 64 +++++++++++++++++ dumux/flux/cctpfa/darcyslaw.hh | 16 ++--- dumux/flux/cctpfa/dispersionflux.hh | 2 +- dumux/flux/cctpfa/fickslaw.hh | 6 +- dumux/flux/cctpfa/fourierslaw.hh | 6 +- .../flux/cctpfa/fourierslawnonequilibrium.hh | 6 +- .../boundary/darcydarcy/couplingmanager.hh | 4 +- .../couplingconditions.hh | 2 +- .../boundary/stokesdarcy/couplingdata.hh | 2 +- .../facet/cellcentered/mpfa/localassembler.hh | 6 +- .../facet/cellcentered/tpfa/darcyslaw.hh | 6 +- .../facet/cellcentered/tpfa/fickslaw.hh | 6 +- .../facet/cellcentered/tpfa/fourierslaw.hh | 6 +- .../fluxvariablescachefiller.hh | 2 +- 16 files changed, 168 insertions(+), 40 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh b/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh index c1bb46d0e3..321f25a9d0 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, @@ -58,7 +59,38 @@ computeMpfaTransmissibility(const IVSubControlVolume& scv, wijk[dir] = vtmv(scvf.unitOuterNormal(), t, scv.nu(dir)); using Extrusion = Extrusion_t; - wijk *= Extrusion::area(scvf)*extrusionFactor;// TODO: No fvGeometry available + wijk *= Extrusion::area(scvf)*extrusionFactor; + wijk /= scv.detX(); + + 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 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; @@ -77,6 +109,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, @@ -88,7 +121,38 @@ computeMpfaTransmissibility(const IVSubControlVolume& scv, wijk[dir] = vtmv(scvf.unitOuterNormal(), t, scv.nu(dir)); using Extrusion = Extrusion_t; - wijk *= Extrusion::area(scvf)*extrusionFactor;// TODO: No fvGeometry available + wijk *= Extrusion::area(scvf)*extrusionFactor; + wijk /= scv.detX(); + + 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 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; diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh index 5b0dc2ea74..0e954a0971 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 83bff63f70..486205b79d 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 finite-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 SubControlVolumeFace, class SubControlVolume, class Tensor > +typename Tensor::field_type computeTpfaTransmissibility(const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf, + const SubControlVolume& scv, + const Tensor& T, + typename SubControlVolume::Traits::Scalar extrusionFactor) +{ + using GlobalPosition = typename 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,37 @@ 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 finite-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 SubControlVolumeFace, + class SubControlVolume, + class Tensor, + typename std::enable_if_t::value, int> = 0 > +Tensor computeTpfaTransmissibility(const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf, + const SubControlVolume &scv, + Tensor t, + typename 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/flux/cctpfa/darcyslaw.hh b/dumux/flux/cctpfa/darcyslaw.hh index 88a96791c2..378e79b674 100644 --- a/dumux/flux/cctpfa/darcyslaw.hh +++ b/dumux/flux/cctpfa/darcyslaw.hh @@ -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,7 +239,7 @@ 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()); @@ -256,8 +256,8 @@ 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? @@ -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,7 +481,7 @@ 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()); @@ -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()); diff --git a/dumux/flux/cctpfa/dispersionflux.hh b/dumux/flux/cctpfa/dispersionflux.hh index e3428afa20..9c0fe1e5d2 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 1a65750bd9..7d681adddb 100644 --- a/dumux/flux/cctpfa/fickslaw.hh +++ b/dumux/flux/cctpfa/fickslaw.hh @@ -241,7 +241,7 @@ 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; @@ -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()); diff --git a/dumux/flux/cctpfa/fourierslaw.hh b/dumux/flux/cctpfa/fourierslaw.hh index 04f2a831f4..d4b3c98f05 100644 --- a/dumux/flux/cctpfa/fourierslaw.hh +++ b/dumux/flux/cctpfa/fourierslaw.hh @@ -183,7 +183,7 @@ public: const auto& insideVolVars = elemVolVars[insideScvIdx]; const auto insideLambda = insideVolVars.effectiveThermalConductivity(); - 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) @@ -201,9 +201,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, 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) diff --git a/dumux/flux/cctpfa/fourierslawnonequilibrium.hh b/dumux/flux/cctpfa/fourierslawnonequilibrium.hh index ddd7fbde5b..d92e1333d2 100644 --- a/dumux/flux/cctpfa/fourierslawnonequilibrium.hh +++ b/dumux/flux/cctpfa/fourierslawnonequilibrium.hh @@ -125,7 +125,7 @@ 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) @@ -140,9 +140,9 @@ 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) diff --git a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh index 09f845989b..2fafb01960 100644 --- a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh +++ b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh @@ -178,8 +178,8 @@ 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(fvGeometry, scvf)*(ti * tj)/(ti + tj); diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh index 6603b58268..74e8145786 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/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh index c768126d90..14221ca9a2 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/facet/cellcentered/mpfa/localassembler.hh b/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh index 7c063771c8..61640013b2 100644 --- a/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh +++ b/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh @@ -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; @@ -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 42c9691ef7..bd7ca6dbfa 100644 --- a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh +++ b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh @@ -264,7 +264,7 @@ class CCTpfaFacetCouplingDarcysLawImpl Date: Wed, 14 Sep 2022 17:07:09 +0200 Subject: [PATCH 14/20] Use extrusion-aware volume calculation for geometries --- dumux/porousmediumflow/2p/griddatatransfer.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/porousmediumflow/2p/griddatatransfer.hh b/dumux/porousmediumflow/2p/griddatatransfer.hh index 43d4d33547..9bd37686c3 100644 --- a/dumux/porousmediumflow/2p/griddatatransfer.hh +++ b/dumux/porousmediumflow/2p/griddatatransfer.hh @@ -241,7 +241,7 @@ public: if (!isBox) elemSol[0] /= adaptedValues.count; - const auto elementVolume = volume(element.geometry()); + const auto elementVolume = volume(element.geometry(), Extrusion{}); for (const auto& scv : scvs(fvGeometry)) { VolumeVariables volVars; -- GitLab From 2cd21a87e8a78516c73f72ac14e614239e2b3daf Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Fri, 16 Sep 2022 10:52:42 +0200 Subject: [PATCH 15/20] Fixup parameters and descriptions --- .../mpfa/computetransmissibility.hh | 2 ++ .../tpfa/computetransmissibility.hh | 22 +++++++++---------- .../navierstokes/staggered/fluxvariables.hh | 1 + 3 files changed, 13 insertions(+), 12 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh b/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh index 321f25a9d0..ed0ddb0a44 100644 --- a/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh +++ b/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh @@ -72,6 +72,7 @@ computeMpfaTransmissibility(const IVSubControlVolume& scv, * 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 @@ -134,6 +135,7 @@ computeMpfaTransmissibility(const IVSubControlVolume& scv, * 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 diff --git a/dumux/discretization/cellcentered/tpfa/computetransmissibility.hh b/dumux/discretization/cellcentered/tpfa/computetransmissibility.hh index 486205b79d..793dea331e 100644 --- a/dumux/discretization/cellcentered/tpfa/computetransmissibility.hh +++ b/dumux/discretization/cellcentered/tpfa/computetransmissibility.hh @@ -68,20 +68,20 @@ typename Tensor::field_type computeTpfaTransmissibility(const SubControlVolumeFa * sub-control volume face stemming from a given sub-control * volume with corresponding tensor T. * - * \param fvGeometry The finite-volume geometry + * \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 SubControlVolumeFace, class SubControlVolume, class Tensor > +template< class FVElementGeometry, class Tensor > typename Tensor::field_type computeTpfaTransmissibility(const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf, - const SubControlVolume& scv, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const typename FVElementGeometry::SubControlVolume& scv, const Tensor& T, - typename SubControlVolume::Traits::Scalar extrusionFactor) + typename FVElementGeometry::SubControlVolume::Traits::Scalar extrusionFactor) { - using GlobalPosition = typename SubControlVolumeFace::Traits::GlobalPosition; + using GlobalPosition = typename FVElementGeometry::SubControlVolumeFace::Traits::GlobalPosition; GlobalPosition Knormal; T.mv(scvf.unitOuterNormal(), Knormal); @@ -128,22 +128,20 @@ Tensor computeTpfaTransmissibility(const SubControlVolumeFace& scvf, * sub-control volume face stemming from a given sub-control * volume for the case where T is just a scalar * - * \param fvGeometry The finite-volume geometry + * \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 SubControlVolumeFace, - class SubControlVolume, class Tensor, typename std::enable_if_t::value, int> = 0 > Tensor computeTpfaTransmissibility(const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf, - const SubControlVolume &scv, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const typename FVElementGeometry::SubControlVolume& scv, Tensor t, - typename SubControlVolumeFace::Traits::Scalar extrusionFactor) + typename FVElementGeometry::SubControlVolumeFace::Traits::Scalar extrusionFactor) { auto distanceVector = scvf.ipGlobal(); distanceVector -= scv.center(); diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh index 0945f55fcd..ce845ed138 100644 --- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh @@ -139,6 +139,7 @@ public: /*! * \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 -- GitLab From 8a2a6385458029fc834b3937b3f1b7986f25a656 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Fri, 16 Sep 2022 10:53:55 +0200 Subject: [PATCH 16/20] [disc] Decouple template arguments and resolve unused variables --- dumux/discretization/extrusion.hh | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/dumux/discretization/extrusion.hh b/dumux/discretization/extrusion.hh index 8e592f3e01..aa30f7ffbd 100644 --- a/dumux/discretization/extrusion.hh +++ b/dumux/discretization/extrusion.hh @@ -49,11 +49,11 @@ struct NoExtrusion { return scv.volume(); } template - static constexpr auto area(const FVGeo& fvGeometry, const SCVF& scvf) + static constexpr auto area(const FVGeo&, const SCVF& scvf) { return scvf.area(); } template - static constexpr auto volume(const FVGeo& fvGeometry, const SCV& scv) + 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!"); -- GitLab From d482db8bb59b1eae84ac56de676669a114e1efcb Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Fri, 16 Sep 2022 10:55:09 +0200 Subject: [PATCH 17/20] [fixup] Unify indentation --- dumux/multidomain/facet/box/subcontrolvolumeface.hh | 8 ++++---- dumux/porousmediumflow/boxdfm/subcontrolvolume.hh | 4 ++-- dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh | 8 ++++---- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/dumux/multidomain/facet/box/subcontrolvolumeface.hh b/dumux/multidomain/facet/box/subcontrolvolumeface.hh index 0c024c883c..c3b14ec4bd 100644 --- a/dumux/multidomain/facet/box/subcontrolvolumeface.hh +++ b/dumux/multidomain/facet/box/subcontrolvolumeface.hh @@ -80,8 +80,8 @@ public: , center_(0.0) , unitOuterNormal_(geometryHelper.normal(corners_, scvIndices)) , area_(Dumux::convexPolytopeVolume( - Dune::GeometryTypes::cube(T::dim-1), - [&](unsigned int i){ return corners_[i]; }) + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; }) ) , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) @@ -110,8 +110,8 @@ public: , center_(0.0) , unitOuterNormal_(intersection.centerUnitOuterNormal()) , area_(Dumux::convexPolytopeVolume( - Dune::GeometryTypes::cube(T::dim-1), - [&](unsigned int i){ return corners_[i]; }) + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; }) ) , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) diff --git a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh index 637932dfe7..fb91e85fc4 100644 --- a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh +++ b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh @@ -103,8 +103,8 @@ public: , corners_(geometryHelper.getScvCorners(scvIdx)) , center_(0.0) , volume_(Dumux::convexPolytopeVolume( - Dune::GeometryTypes::cube(T::dim), - [&](unsigned int i){ return corners_[i]; }) + Dune::GeometryTypes::cube(T::dim), + [&](unsigned int i){ return corners_[i]; }) ) , elementIndex_(elementIndex) , vIdxLocal_(scvIdx) diff --git a/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh index 8cc86bf696..9a1a24a9fa 100644 --- a/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh +++ b/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh @@ -105,8 +105,8 @@ public: , center_(0.0) , unitOuterNormal_(geometryHelper.normal(corners_, scvIndices)) , area_(Dumux::convexPolytopeVolume( - Dune::GeometryTypes::cube(T::dim-1), - [&](unsigned int i){ return corners_[i]; }) + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; }) ) , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) @@ -132,8 +132,8 @@ public: , center_(0.0) , unitOuterNormal_(intersection.centerUnitOuterNormal()) , area_(Dumux::convexPolytopeVolume( - Dune::GeometryTypes::cube(T::dim-1), - [&](unsigned int i){ return corners_[i]; }) + Dune::GeometryTypes::cube(T::dim-1), + [&](unsigned int i){ return corners_[i]; }) ) , scvfIndex_(scvfIndex) , scvIndices_(std::move(scvIndices)) -- GitLab From cbba0eebe32b9c5493642ed810d92d613f24a764 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Fri, 16 Sep 2022 10:55:48 +0200 Subject: [PATCH 18/20] [pmflow] Add missing extrusions for adaptive --- dumux/porousmediumflow/2p/griddatatransfer.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/porousmediumflow/2p/griddatatransfer.hh b/dumux/porousmediumflow/2p/griddatatransfer.hh index 9bd37686c3..2b35f708af 100644 --- a/dumux/porousmediumflow/2p/griddatatransfer.hh +++ b/dumux/porousmediumflow/2p/griddatatransfer.hh @@ -329,7 +329,7 @@ public: else if (formulation == p1s0) massCoeffSon = Extrusion::volume(fvGeometry, scv) * volVars.density(phase0Idx) * volVars.porosity(); sol_[scv.dofIndex()][saturationIdx] = - ( Extrusion::volume(fvGeometry, scv)/volume(fatherElement.geometry())*massFather )/massCoeffSon; + ( Extrusion::volume(fvGeometry, scv)/volume(fatherElement.geometry(), Extrusion{})*massFather )/massCoeffSon; } } else @@ -348,7 +348,7 @@ public: scv.dofPosition()); // compute mass & mass coefficients for the scvs (saturations are recalculated at the end) - const auto fatherElementVolume = volume(fatherGeometry); + const auto fatherElementVolume = volume(fatherGeometry, Extrusion{}); for (const auto& scv : scvs(fvGeometry)) { VolumeVariables volVars; -- GitLab From a9e2d3ed879b800142fd16ce8956a7d667752d2c Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Fri, 16 Sep 2022 10:57:22 +0200 Subject: [PATCH 19/20] [ff][staggered] remove unnecessary function parameter --- dumux/freeflow/navierstokes/staggered/fluxvariables.hh | 2 +- dumux/freeflow/navierstokes/staggered/localresidual.hh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh index ce845ed138..9f04ee91af 100644 --- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh @@ -515,13 +515,13 @@ public: * \endverbatim */ FacePrimaryVariables inflowOutflowBoundaryFlux(const Problem& problem, - const Element& element, 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. diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh index 04f85e3008..62adae2bd7 100644 --- a/dumux/freeflow/navierstokes/staggered/localresidual.hh +++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh @@ -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, fvGeometry, elemVolVars, elemFaceVars); + result += fluxVars.inflowOutflowBoundaryFlux(problem, scvf, fvGeometry, elemVolVars, elemFaceVars); } } return result; -- GitLab From 80fddd80c901e897ca24a63a0cff856fdca03266 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Fri, 16 Sep 2022 12:16:47 +0200 Subject: [PATCH 20/20] Add changelog entries for deprecations --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index bd94e7c656..2499c7d10a 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) -- GitLab