diff --git a/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.C b/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.C
index e80016d9360a257e4ead4ccc788314140722d90e..aeca2e05161f8d0fe3d2f0668ea26baaccdafb10 100644
--- a/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.C
+++ b/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.C
@@ -82,10 +82,10 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
if (mesh().moving())
{
- scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
+ scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
- scalarField SS0 = mesh().S() + mesh().S0();
- scalarField S0S00 = mesh().S0() + mesh().S00();
+ scalarField SS0(mesh().S() + mesh().S0());
+ scalarField S0S00(mesh().S0() + mesh().S00());
tmp<GeometricField<Type, faPatchField, areaMesh>> tdt2dt2
(
@@ -148,10 +148,10 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
if (mesh().moving())
{
- scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
+ scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
- scalarField SS0 = mesh().S() + mesh().S0();
- scalarField S0S00 = mesh().S0() + mesh().S00();
+ scalarField SS0(mesh().S() + mesh().S0());
+ scalarField S0S00(mesh().S0() + mesh().S00());
return tmp<GeometricField<Type, faPatchField, areaMesh>>
(
@@ -228,13 +228,15 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
{
scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
- scalarField SS0rhoRho0 =
- (mesh().S() + mesh().S0())
- *rho.value();
+ scalarField SS0rhoRho0
+ (
+ (mesh().S() + mesh().S0())*rho.value()
+ );
- scalarField S0S00rho0Rho00 =
- (mesh().S0() + mesh().S00())
- *rho.value();
+ scalarField S0S00rho0Rho00
+ (
+ (mesh().S0() + mesh().S00())*rho.value()
+ );
return tmp<GeometricField<Type, faPatchField, areaMesh>>
(
@@ -489,7 +491,6 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
scalar halfRdeltaT2 = 0.5*rDeltaT2;
scalarField SS0(mesh().S() + mesh().S0());
-
scalarField S0S00(mesh().S0() + mesh().S00());
fam.diag() = rho.value()*(coefft*halfRdeltaT2)*SS0;
diff --git a/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.H b/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.H
index 97d501734b5350f8da1ff9550e92885ab3380a04..f778368530a17784fd1835fd6bfeb8fe35fc3b8e 100644
--- a/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.H
+++ b/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.H
@@ -76,10 +76,8 @@ class boundedBackwardFaDdtScheme
{
return GREAT;
}
- else
- {
- return deltaT0_();
- }
+
+ return deltaT0_();
}
diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
index 656187479e7bb822869d50d053f3142e97912ec0..118c60e9f2cd367e3b493a31ba36e86e5c8d3345 100644
--- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
+++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
@@ -91,10 +91,10 @@ leastSquaresFaGrad<Type>::calcGrad
forAll(own, edgei)
{
- label ownEdgeI = own[edgei];
- label neiEdgeI = nei[edgei];
+ const label ownEdgeI = own[edgei];
+ const label neiEdgeI = nei[edgei];
- Type deltaVsf = vsf[neiEdgeI] - vsf[ownEdgeI];
+ const Type deltaVsf(vsf[neiEdgeI] - vsf[ownEdgeI]);
lsGrad[ownEdgeI] += ownLs[edgei]*deltaVsf;
lsGrad[neiEdgeI] -= neiLs[edgei]*deltaVsf;
@@ -103,35 +103,23 @@ leastSquaresFaGrad<Type>::calcGrad
// Boundary edges
forAll(vsf.boundaryField(), patchi)
{
- const faePatchVectorField& patchOwnLs = ownLs.boundaryField()[patchi];
+ const faPatchField<Type>& bf = vsf.boundaryField()[patchi];
+ const Field<Type>& vsfp =
+ (
+ bf.coupled()
+ ? bf.patchNeighbourField().cref()
+ : const_cast<faPatchField<Type>&>(bf)
+ );
+
+ const faePatchVectorField& ownLsp = ownLs.boundaryField()[patchi];
const labelUList& edgeFaces =
lsGrad.boundaryField()[patchi].patch().edgeFaces();
- if (vsf.boundaryField()[patchi].coupled())
- {
- Field<Type> neiVsf
- (
- vsf.boundaryField()[patchi].patchNeighbourField()
- );
-
- forAll(neiVsf, patchEdgeI)
- {
- lsGrad[edgeFaces[patchEdgeI]] +=
- patchOwnLs[patchEdgeI]
- *(neiVsf[patchEdgeI] - vsf[edgeFaces[patchEdgeI]]);
- }
- }
- else
+ forAll(vsfp, pEdgei)
{
- const faPatchField<Type>& patchVsf = vsf.boundaryField()[patchi];
-
- forAll(patchVsf, patchEdgeI)
- {
- lsGrad[edgeFaces[patchEdgeI]] +=
- patchOwnLs[patchEdgeI]
- *(patchVsf[patchEdgeI] - vsf[edgeFaces[patchEdgeI]]);
- }
+ lsGrad[edgeFaces[pEdgei]] +=
+ ownLsp[pEdgei]*(vsfp[pEdgei] - vsf[edgeFaces[pEdgei]]);
}
}
diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
index 3aa9c4d205d9a36eab6e8eb38fb8c90f5f744a5e..a68d4f63ec717ad2578af29e386c879b88dc9d5c 100644
--- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
+++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
@@ -91,22 +91,22 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::calcGrad
const areaVectorField& C = mesh.areaCentres();
const edgeVectorField& Cf = mesh.edgeCentres();
- // create limiter
+ // Create limiter field
scalarField limiter(vsf.internalField().size(), 1.0);
- scalar rk = (1.0/k_ - 1.0);
+ const scalar rk = (1.0/k_ - 1.0);
forAll(owner, edgei)
{
- label own = owner[edgei];
- label nei = neighbour[edgei];
+ const label own = owner[edgei];
+ const label nei = neighbour[edgei];
- scalar vsfOwn = vsf[own];
- scalar vsfNei = vsf[nei];
+ const scalar vsfOwn = vsf[own];
+ const scalar vsfNei = vsf[nei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
- scalar maxMinEdge = rk*(maxEdge - minEdge);
+ const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
@@ -114,7 +114,8 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::calcGrad
limitEdge
(
limiter[own],
- maxEdge - vsfOwn, minEdge - vsfOwn,
+ maxEdge - vsfOwn,
+ minEdge - vsfOwn,
(Cf[edgei] - C[own]) & g[own]
);
@@ -122,67 +123,53 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::calcGrad
limitEdge
(
limiter[nei],
- maxEdge - vsfNei, minEdge - vsfNei,
+ maxEdge - vsfNei,
+ minEdge - vsfNei,
(Cf[edgei] - C[nei]) & g[nei]
);
}
- const areaScalarField::Boundary& bsf = vsf.boundaryField();
+ // Lambda expression to update limiter for boundary edges
+ auto updateLimiter = [&](const label patchi, const scalarField& fld) -> void
+ {
+ const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
+ const vectorField& pCf = Cf.boundaryField()[patchi];
+
+ forAll(pOwner, edgei)
+ {
+ const label own = pOwner[edgei];
+
+ const scalar vsfOwn = vsf[own];
+ const scalar vsfNei = fld[edgei];
+
+ scalar maxEdge = max(vsfOwn, vsfNei);
+ scalar minEdge = min(vsfOwn, vsfNei);
+ const scalar maxMinEdge = rk*(maxEdge - minEdge);
+ maxEdge += maxMinEdge;
+ minEdge -= maxMinEdge;
+
+ limitEdge
+ (
+ limiter[own],
+ maxEdge - vsfOwn,
+ minEdge - vsfOwn,
+ (pCf[edgei] - C[own]) & g[own]
+ );
+ }
+ };
+ const areaScalarField::Boundary& bsf = vsf.boundaryField();
forAll(bsf, patchi)
{
const faPatchScalarField& psf = bsf[patchi];
- const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
- const vectorField& pCf = Cf.boundaryField()[patchi];
-
if (psf.coupled())
{
- const scalarField psfNei(psf.patchNeighbourField());
-
- forAll(pOwner, pEdgei)
- {
- label own = pOwner[pEdgei];
-
- scalar vsfOwn = vsf[own];
- scalar vsfNei = psfNei[pEdgei];
-
- scalar maxEdge = max(vsfOwn, vsfNei);
- scalar minEdge = min(vsfOwn, vsfNei);
- scalar maxMinEdge = rk*(maxEdge - minEdge);
- maxEdge += maxMinEdge;
- minEdge -= maxMinEdge;
-
- limitEdge
- (
- limiter[own],
- maxEdge - vsfOwn, minEdge - vsfOwn,
- (pCf[pEdgei] - C[own]) & g[own]
- );
- }
+ updateLimiter(patchi, psf.patchNeighbourField());
}
else if (psf.fixesValue())
{
- forAll(pOwner, pEdgei)
- {
- label own = pOwner[pEdgei];
-
- scalar vsfOwn = vsf[own];
- scalar vsfNei = psf[pEdgei];
-
- scalar maxEdge = max(vsfOwn, vsfNei);
- scalar minEdge = min(vsfOwn, vsfNei);
- scalar maxMinEdge = rk*(maxEdge - minEdge);
- maxEdge += maxMinEdge;
- minEdge -= maxMinEdge;
-
- limitEdge
- (
- limiter[own],
- maxEdge - vsfOwn, minEdge - vsfOwn,
- (pCf[pEdgei] - C[own]) & g[own]
- );
- }
+ updateLimiter(patchi, psf);
}
}
@@ -226,35 +213,33 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::calcGrad
const areaVectorField& C = mesh.areaCentres();
const edgeVectorField& Cf = mesh.edgeCentres();
- // create limiter
+ // Create limiter
scalarField limiter(vvf.internalField().size(), 1.0);
- scalar rk = (1.0/k_ - 1.0);
+ const scalar rk = (1.0/k_ - 1.0);
forAll(owner, edgei)
{
- label own = owner[edgei];
- label nei = neighbour[edgei];
-
- vector vvfOwn = vvf[own];
- vector vvfNei = vvf[nei];
+ const label own = owner[edgei];
+ const label nei = neighbour[edgei];
// owner side
- vector gradf = (Cf[edgei] - C[own]) & g[own];
+ vector gradf((Cf[edgei] - C[own]) & g[own]);
- scalar vsfOwn = gradf & vvfOwn;
- scalar vsfNei = gradf & vvfNei;
+ scalar vsfOwn = gradf & vvf[own];
+ scalar vsfNei = gradf & vvf[nei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
- scalar maxMinEdge = rk*(maxEdge - minEdge);
+ const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
limitEdge
(
limiter[own],
- maxEdge - vsfOwn, minEdge - vsfOwn,
+ maxEdge - vsfOwn,
+ minEdge - vsfOwn,
magSqr(gradf)
);
@@ -262,8 +247,8 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::calcGrad
// neighbour side
gradf = (Cf[edgei] - C[nei]) & g[nei];
- vsfOwn = gradf & vvfOwn;
- vsfNei = gradf & vvfNei;
+ vsfOwn = gradf & vvf[own];
+ vsfNei = gradf & vvf[nei];
maxEdge = max(vsfOwn, vsfNei);
minEdge = min(vsfOwn, vsfNei);
@@ -271,78 +256,57 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::calcGrad
limitEdge
(
limiter[nei],
- maxEdge - vsfNei, minEdge - vsfNei,
+ maxEdge - vsfNei,
+ minEdge - vsfNei,
magSqr(gradf)
);
}
- const areaVectorField::Boundary& bvf = vvf.boundaryField();
-
- forAll(bvf, patchi)
+ // Lambda expression to update limiter for boundary edges
+ auto updateLimiter = [&](const label patchi, const vectorField& fld) -> void
{
- const faPatchVectorField& psf = bvf[patchi];
-
const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
const vectorField& pCf = Cf.boundaryField()[patchi];
- if (psf.coupled())
+ forAll(pOwner, edgei)
{
- const vectorField psfNei(psf.patchNeighbourField());
-
- forAll(pOwner, pEdgei)
- {
- label own = pOwner[pEdgei];
-
- vector vvfOwn = vvf[own];
- vector vvfNei = psfNei[pEdgei];
-
- vector gradf = (pCf[pEdgei] - C[own]) & g[own];
+ const label own = pOwner[edgei];
+
+ const vector gradf((pCf[edgei] - C[own]) & g[own]);
+
+ const scalar vsfOwn = gradf & vvf[own];
+ const scalar vsfNei = gradf & fld[edgei];
+
+ scalar maxEdge = max(vsfOwn, vsfNei);
+ scalar minEdge = min(vsfOwn, vsfNei);
+ const scalar maxMinEdge = rk*(maxEdge - minEdge);
+ maxEdge += maxMinEdge;
+ minEdge -= maxMinEdge;
+
+ limitEdge
+ (
+ limiter[own],
+ maxEdge - vsfOwn,
+ minEdge - vsfOwn,
+ magSqr(gradf)
+ );
+ }
+ };
- scalar vsfOwn = gradf & vvfOwn;
- scalar vsfNei = gradf & vvfNei;
- scalar maxEdge = max(vsfOwn, vsfNei);
- scalar minEdge = min(vsfOwn, vsfNei);
- scalar maxMinEdge = rk*(maxEdge - minEdge);
- maxEdge += maxMinEdge;
- minEdge -= maxMinEdge;
+ const areaVectorField::Boundary& bvf = vvf.boundaryField();
+ forAll(bvf, patchi)
+ {
+ const faPatchVectorField& psf = bvf[patchi];
- limitEdge
- (
- limiter[own],
- maxEdge - vsfOwn, minEdge - vsfOwn,
- magSqr(gradf)
- );
- }
+ if (psf.coupled())
+ {
+ updateLimiter(patchi, psf.patchNeighbourField());
}
else if (psf.fixesValue())
{
- forAll(pOwner, pEdgei)
- {
- label own = pOwner[pEdgei];
-
- vector vvfOwn = vvf[own];
- vector vvfNei = psf[pEdgei];
-
- vector gradf = (pCf[pEdgei] - C[own]) & g[own];
-
- scalar vsfOwn = gradf & vvfOwn;
- scalar vsfNei = gradf & vvfNei;
-
- scalar maxEdge = max(vsfOwn, vsfNei);
- scalar minEdge = min(vsfOwn, vsfNei);
- scalar maxMinEdge = rk*(maxEdge - minEdge);
- maxEdge += maxMinEdge;
- minEdge -= maxMinEdge;
-
- limitEdge
- (
- limiter[own],
- maxEdge - vsfOwn, minEdge - vsfOwn,
- magSqr(gradf)
- );
- }
+ updateLimiter(patchi, psf);
}
}
diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
index da2fb1bf7b075c7a2a40d1492300647461aab67b..eb0cc9bcbd05eaf9f37004519f3e0f0d54db7ef6 100644
--- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
+++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
@@ -118,11 +118,11 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::calcGrad
forAll(owner, facei)
{
- label own = owner[facei];
- label nei = neighbour[facei];
+ const label own = owner[facei];
+ const label nei = neighbour[facei];
- scalar vsfOwn = vsf[own];
- scalar vsfNei = vsf[nei];
+ const scalar vsfOwn = vsf[own];
+ const scalar vsfNei = vsf[nei];
maxVsf[own] = max(maxVsf[own], vsfNei);
minVsf[own] = min(minVsf[own], vsfNei);
@@ -132,37 +132,33 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::calcGrad
}
- const areaScalarField::Boundary& bsf = vsf.boundaryField();
+ // Lambda expression to update limiter for boundary edges
+ auto updateLimiter = [&](const label patchi, const scalarField& fld) -> void
+ {
+ const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
+
+ forAll(pOwner, facei)
+ {
+ const label own = pOwner[facei];
+ const scalar vsf = fld[facei];
+
+ maxVsf[own] = max(maxVsf[own], vsf);
+ minVsf[own] = min(minVsf[own], vsf);
+ }
+ };
+ const areaScalarField::Boundary& bsf = vsf.boundaryField();
forAll(bsf, patchi)
{
const faPatchScalarField& psf = bsf[patchi];
- const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
-
if (psf.coupled())
{
- scalarField psfNei(psf.patchNeighbourField());
-
- forAll(pOwner, pFacei)
- {
- label own = pOwner[pFacei];
- scalar vsfNei = psfNei[pFacei];
-
- maxVsf[own] = max(maxVsf[own], vsfNei);
- minVsf[own] = min(minVsf[own], vsfNei);
- }
+ updateLimiter(patchi, psf.patchNeighbourField());
}
else
{
- forAll(pOwner, pFacei)
- {
- label own = pOwner[pFacei];
- scalar vsfNei = psf[pFacei];
-
- maxVsf[own] = max(maxVsf[own], vsfNei);
- minVsf[own] = min(minVsf[own], vsfNei);
- }
+ updateLimiter(patchi, psf);
}
}
@@ -171,19 +167,19 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::calcGrad
if (k_ < 1.0)
{
- scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
+ const scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
maxVsf += maxMinVsf;
minVsf -= maxMinVsf;
}
- // create limiter
+ // Create limiter
scalarField limiter(vsf.internalField().size(), 1.0);
forAll(owner, facei)
{
- label own = owner[facei];
- label nei = neighbour[facei];
+ const label own = owner[facei];
+ const label nei = neighbour[facei];
// owner side
limitEdge
@@ -211,7 +207,7 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::calcGrad
forAll(pOwner, pFacei)
{
- label own = pOwner[pFacei];
+ const label own = pOwner[pFacei];
limitEdge
(
@@ -268,8 +264,8 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
forAll(owner, facei)
{
- label own = owner[facei];
- label nei = neighbour[facei];
+ const label own = owner[facei];
+ const label nei = neighbour[facei];
const vector& vsfOwn = vsf[own];
const vector& vsfNei = vsf[nei];
@@ -282,36 +278,33 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
}
- const areaVectorField::Boundary& bsf = vsf.boundaryField();
+ // Lambda expression to update limiter for boundary edges
+ auto updateLimiter = [&](const label patchi, const vectorField& fld) -> void
+ {
+ const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
+ forAll(pOwner, facei)
+ {
+ const label own = pOwner[facei];
+ const vector& vsf = fld[facei];
+
+ maxVsf[own] = max(maxVsf[own], vsf);
+ minVsf[own] = min(minVsf[own], vsf);
+ }
+ };
+
+ const areaVectorField::Boundary& bsf = vsf.boundaryField();
forAll(bsf, patchi)
{
const faPatchVectorField& psf = bsf[patchi];
- const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
if (psf.coupled())
{
- vectorField psfNei(psf.patchNeighbourField());
-
- forAll(pOwner, pFacei)
- {
- label own = pOwner[pFacei];
- const vector& vsfNei = psfNei[pFacei];
-
- maxVsf[own] = max(maxVsf[own], vsfNei);
- minVsf[own] = min(minVsf[own], vsfNei);
- }
+ updateLimiter(patchi, psf.patchNeighbourField());
}
else
{
- forAll(pOwner, pFacei)
- {
- label own = pOwner[pFacei];
- const vector& vsfNei = psf[pFacei];
-
- maxVsf[own] = max(maxVsf[own], vsfNei);
- minVsf[own] = min(minVsf[own], vsfNei);
- }
+ updateLimiter(patchi, psf);
}
}
@@ -320,7 +313,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
if (k_ < 1.0)
{
- vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
+ const vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
maxVsf += maxMinVsf;
minVsf -= maxMinVsf;
@@ -329,13 +322,13 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
}
- // create limiter
+ // Create limiter
vectorField limiter(vsf.internalField().size(), vector::one);
forAll(owner, facei)
{
- label own = owner[facei];
- label nei = neighbour[facei];
+ const label own = owner[facei];
+ const label nei = neighbour[facei];
// owner side
limitEdge
@@ -363,7 +356,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
forAll(pOwner, pFacei)
{
- label own = pOwner[pFacei];
+ const label own = pOwner[pFacei];
limitEdge
(
diff --git a/src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H b/src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H
index 18b4973828f7747b7929c2c2b21079b1a5730df1..6b3eb5f3be418dd0054f5aac358063c7f4dfd4a3 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H
+++ b/src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H
@@ -76,7 +76,7 @@ public:
// Rescale k_ to be >= 0 and <= 0.5 (TVD conformant)
// and avoid the /0 when k_ = 0
- k_ = max(k_/2.0, SMALL);
+ k_ = max(0.5*k_, SMALL);
}
@@ -91,30 +91,19 @@ public:
const vector& d
) const
{
- scalar magd = mag(d);
- vector dHat = d/mag(d);
+ const vector dHat(normalised(d));
- scalar gradf = (phiN - phiP)/magd;
-
- scalar gradcf;
- scalar udWeight;
-
- if (faceFlux > 0)
- {
- gradcf = dHat & gradcP;
- udWeight = 1;
- }
- else
- {
- gradcf = dHat & gradcN;
- udWeight = 0;
- }
+ // Choose gradc based on faceFlux
+ const vector& gradcPN = (faceFlux > 0) ? gradcP : gradcN;
+ const scalar udWeight = (faceFlux > 0) ? 1 : 0;
// Stabilise for division
- gradcf = stabilise(gradcf, SMALL);
+ const scalar gradcf = stabilise(dHat & gradcPN, SMALL);
+
+ const scalar gradf = (phiN - phiP)/mag(d);
- scalar phict = 1 - 0.5*gradf/gradcf;
- scalar limiter = clamp(phict/k_, zero_one{});
+ const scalar phict = 1 - 0.5*gradf/gradcf;
+ const scalar limiter = clamp(phict/k_, zero_one{});
return lerp(udWeight, cdWeight, limiter);
}