From 9aff74aaafdaa15454ff629ec00e7c2d95642b08 Mon Sep 17 00:00:00 2001 From: Andrew Heather <> Date: Wed, 20 Sep 2017 13:55:42 +0100 Subject: [PATCH] COMP: Updated to compile with Clang 3.7.1 --- .../liquidFilmFoam/createFaFields.H | 4 +- .../liquidFilmFoam/liquidFilmFoam.C | 2 +- .../liquidFilmFoam/surfaceCourantNo.H | 7 +- src/finiteArea/faMatrices/faMatrix/faMatrix.C | 10 +- .../faMatrices/faMatrix/faMatrixSolve.C | 14 +- .../faScalarMatrix/faScalarMatrix.C | 12 +- src/finiteArea/faMesh/faMesh.C | 8 +- src/finiteArea/faMesh/faMesh.H | 9 +- .../faMesh/faMeshDemandDrivenData.C | 11 +- .../faPatches/basic/coupled/coupledFaPatch.H | 2 +- .../constraint/cyclic/cyclicFaPatch.C | 8 +- .../faMesh/faPatches/faPatch/faPatch.C | 2 +- .../basicSymmetry/basicSymmetryFaPatchField.C | 7 +- .../constraint/wedge/wedgeFaPatchField.C | 2 +- .../gaussFaConvectionScheme.C | 6 +- .../boundedBackwardFaDdtScheme.C | 120 ++++++++++-------- .../gradSchemes/gaussFaGrad/gaussFaGrad.C | 6 +- .../leastSquaresFaGrad/leastSquaresFaGrad.C | 10 +- .../leastSquaresFaVectors.C | 8 +- .../leastSquaresFaVectors.H | 2 +- .../edgeLimitedFaGrad/edgeLimitedFaGrads.C | 4 +- .../faceLimitedFaGrad/faceLimitedFaGrads.C | 10 +- .../gaussFaLaplacianScheme.C | 2 +- .../lnGradSchemes/fourthLnGrad/fourthLnGrad.C | 2 +- .../limitedLnGrad/limitedLnGrad.C | 8 +- .../edgeInterpolation/edgeInterpolation.C | 17 ++- .../edgeInterpolationScheme.C | 10 +- .../schemes/NVDscheme/faNVDscheme.C | 62 +++++---- 28 files changed, 207 insertions(+), 158 deletions(-) diff --git a/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H b/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H index dba96b7fdc..bb5ef10419 100644 --- a/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H +++ b/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H @@ -57,8 +57,8 @@ const areaVectorField& Ns = aMesh.faceAreaNormals(); - areaVectorField Gs = g - Ns*(Ns & g); - areaScalarField Gn = mag(g - Gs); + areaVectorField Gs(g - Ns*(Ns & g)); + areaScalarField Gn(mag(g - Gs)); // Mass source areaScalarField Sm diff --git a/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C b/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C index ce1b06fb8e..192e69872f 100644 --- a/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C +++ b/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C @@ -86,7 +86,7 @@ int main(int argc, char *argv[]) UsEqn.relax(); solve(UsEqn == - fac::grad(ps*h)/rhol + ps*fac::grad(h)/rhol); - areaScalarField UsA = UsEqn.A(); + areaScalarField UsA(UsEqn.A()); Us = UsEqn.H()/UsA; Us.correctBoundaryConditions(); diff --git a/applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H b/applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H index b73cadb0b2..9e58e23282 100644 --- a/applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H +++ b/applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H @@ -41,8 +41,10 @@ scalar velMag = 0.0; if (aMesh.nInternalEdges()) { - edgeScalarField SfUfbyDelta = - aMesh.edgeInterpolation::deltaCoeffs()*mag(phis); + edgeScalarField SfUfbyDelta + ( + aMesh.edgeInterpolation::deltaCoeffs()*mag(phis) + ); CoNum = max(SfUfbyDelta/aMesh.magLe()) .value()*runTime.deltaT().value(); @@ -57,4 +59,5 @@ Info<< "Courant Number mean: " << meanCoNum << " max: " << CoNum << " velocity magnitude: " << velMag << endl; + // ************************************************************************* // diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.C b/src/finiteArea/faMatrices/faMatrix/faMatrix.C index a1323bb109..208a1b72ed 100644 --- a/src/finiteArea/faMatrices/faMatrix/faMatrix.C +++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.C @@ -238,7 +238,7 @@ faMatrix::faMatrix // Update the boundary coefficients of psi without changing its event No. GeometricField& psiRef = - const_cast&>(psi_); + const_cast&>(psi_); label currentStatePsi = psiRef.eventNo(); psiRef.boundaryFieldRef().updateCoeffs(); @@ -627,7 +627,7 @@ tmp > faMatrix::H() const // Loop over field components for (direction cmpt=0; cmpt::operator=(const faMatrix& famv) { faceFluxCorrectionPtr_ = new GeometricField - (*famv.faceFluxCorrectionPtr_); + (*famv.faceFluxCorrectionPtr_); } } @@ -841,7 +841,7 @@ void faMatrix::operator-=(const faMatrix& famv) { faceFluxCorrectionPtr_ = new GeometricField - (-*famv.faceFluxCorrectionPtr_); + (-*famv.faceFluxCorrectionPtr_); } } @@ -930,7 +930,7 @@ void faMatrix::operator*= forAll(boundaryCoeffs_, patchI) { - scalarField psf = vsf.boundaryField()[patchI].patchInternalField(); + const scalarField psf(vsf.boundaryField()[patchI].patchInternalField()); internalCoeffs_[patchI] *= psf; boundaryCoeffs_[patchI] *= psf; } diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C b/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C index 701fb8de49..c5f6780a63 100644 --- a/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C +++ b/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C @@ -68,9 +68,9 @@ SolverPerformance faMatrix::solve(const dictionary& solverControls) psi_.name() ); - scalarField saveDiag = diag(); + scalarField saveDiag(diag()); - Field source = source_; + Field source(source_); addBoundarySource(source); // Make a copy of interfaces: no longer a reference @@ -80,16 +80,16 @@ SolverPerformance faMatrix::solve(const dictionary& solverControls) // Cast into a non-const to solve. HJ, 6/May/2016 GeometricField& psi = - const_cast&>(psi_); + const_cast&>(psi_); for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++) { // copy field and source - scalarField psiCmpt = psi_.primitiveField().component(cmpt); + scalarField psiCmpt(psi_.primitiveField().component(cmpt)); addBoundaryDiag(diag(), cmpt); - scalarField sourceCmpt = source.component(cmpt); + scalarField sourceCmpt(source.component(cmpt)); FieldField bouCoeffsCmpt ( @@ -186,7 +186,7 @@ template tmp > faMatrix::residual() const { tmp > tres(source_); - Field& res = tres(); + Field& res = tres().ref(); addBoundarySource(res); @@ -198,7 +198,7 @@ tmp > faMatrix::residual() const // Loop over field components for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++) { - scalarField psiCmpt = psi_.internalField().component(cmpt); + scalarField psiCmpt(psi_.internalField().component(cmpt)); scalarField boundaryDiagCmpt(psi_.size(), 0.0); addBoundaryDiag(boundaryDiagCmpt, cmpt); diff --git a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C index f7e8a81e75..f1cacedd58 100644 --- a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C +++ b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C @@ -49,11 +49,9 @@ void faMatrix::setComponentReference const scalar value ) { - const labelUList& faceLabels = - psi_.mesh().boundary()[patchI].edgeFaces(); + const labelUList& faceLabels = psi_.mesh().boundary()[patchI].edgeFaces(); - internalCoeffs_[patchI][edgeI] += - diag()[faceLabels[edgeI]]; + internalCoeffs_[patchI][edgeI] += diag()[faceLabels[edgeI]]; boundaryCoeffs_[patchI][edgeI] = value; } @@ -75,10 +73,10 @@ solverPerformance faMatrix::solve GeometricField& psi = const_cast&>(psi_); - scalarField saveDiag = diag(); + scalarField saveDiag(diag()); addBoundaryDiag(diag(), 0); - scalarField totalSource = source_; + scalarField totalSource(source_); addBoundarySource(totalSource, 0); // Solver call @@ -148,7 +146,7 @@ tmp faMatrix::H() const zeroGradientFaPatchScalarField::typeName ) ); - areaScalarField Hphi = tHphi(); + areaScalarField& Hphi = tHphi.ref(); Hphi.primitiveFieldRef() = (lduMatrix::H(psi_.primitiveField()) + source_); addBoundarySource(Hphi.primitiveFieldRef()); diff --git a/src/finiteArea/faMesh/faMesh.C b/src/finiteArea/faMesh/faMesh.C index d103934c92..a14fc8e417 100644 --- a/src/finiteArea/faMesh/faMesh.C +++ b/src/finiteArea/faMesh/faMesh.C @@ -188,7 +188,7 @@ void Foam::faMesh::clearOut() const Foam::faMesh::faMesh(const polyMesh& pMesh) : GeoMesh(pMesh), - MeshObject(pMesh), + MeshObject(pMesh), edgeInterpolation(*this), faceLabels_ ( @@ -282,7 +282,7 @@ Foam::faMesh::faMesh ) : GeoMesh(pMesh), - MeshObject(pMesh), + MeshObject(pMesh), edgeInterpolation(*this), faceLabels_ ( @@ -347,7 +347,7 @@ Foam::faMesh::faMesh ) : GeoMesh(pMesh), - MeshObject(pMesh), + MeshObject(pMesh), edgeInterpolation(*this), faceLabels_ ( @@ -831,7 +831,7 @@ Foam::faMesh::faMesh ) : GeoMesh(pMesh), - MeshObject(pMesh), + MeshObject(pMesh), edgeInterpolation(*this), faceLabels_ ( diff --git a/src/finiteArea/faMesh/faMesh.H b/src/finiteArea/faMesh/faMesh.H index e54390e689..7e7ec0b16b 100644 --- a/src/finiteArea/faMesh/faMesh.H +++ b/src/finiteArea/faMesh/faMesh.H @@ -76,7 +76,7 @@ class faMeshMapper; class faMesh : public GeoMesh, - public MeshObject, + public MeshObject, public lduMesh, public edgeInterpolation { @@ -317,7 +317,12 @@ public: const polyMesh& mesh() const { return - MeshObject::mesh(); + MeshObject + < + polyMesh, + Foam::UpdateableMeshObject, + faMesh + >::mesh(); } //- Return the local mesh directory (dbDir()/meshSubDir) diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C index a51a94b338..d9cd4b7811 100644 --- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C +++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C @@ -683,8 +683,7 @@ void faMesh::calcFaceCurvatures() const // fac::edgeIntegrate(Le()*edgeLengthCorrection()) // &faceAreaNormals(); - areaVectorField kN = - fac::edgeIntegrate(Le()*edgeLengthCorrection()); + areaVectorField kN(fac::edgeIntegrate(Le()*edgeLengthCorrection())); faceCurvatures = sign(kN&faceAreaNormals())*mag(kN); } @@ -795,11 +794,9 @@ void faMesh::calcEdgeTransformTensors() const const labelUList& edgeFaces = boundary()[patchI].edgeFaces(); - vectorField ngbCf = - Cf.boundaryField()[patchI].patchNeighbourField(); + vectorField ngbCf(Cf.boundaryField()[patchI].patchNeighbourField()); - vectorField ngbNf = - Nf.boundaryField()[patchI].patchNeighbourField(); + vectorField ngbNf(Nf.boundaryField()[patchI].patchNeighbourField()); forAll(edgeFaces, edgeI) { @@ -1304,7 +1301,7 @@ void faMesh::calcPointAreaNormals() const labelList patchPoints = boundary()[patchI].pointLabels(); - vectorField N = boundary()[patchI].ngbPolyPatchPointNormals(); + const vectorField N(boundary()[patchI].ngbPolyPatchPointNormals()); forAll (patchPoints, pointI) { diff --git a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H index 43bd4198e3..323e4df64b 100644 --- a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H +++ b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H @@ -255,7 +255,7 @@ public: virtual void initInternalFieldTransfer ( const Pstream::commsTypes commsType, - labelUList& iF + const labelUList& iF ) const {} diff --git a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C index 5156c97e2d..1876c2e61c 100644 --- a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C +++ b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C @@ -61,7 +61,7 @@ void Foam::cyclicFaPatch::calcTransforms() vectorField half0Normals(size()/2); vectorField half1Normals(size()/2); - vectorField eN = edgeNormals()*magEdgeLengths(); + const vectorField eN(edgeNormals()*magEdgeLengths()); scalar maxMatchError = 0; label errorEdge = -1; @@ -159,7 +159,7 @@ void cyclicFaPatch::makeWeights(scalarField& w) const { const scalarField& magL = magEdgeLengths(); - scalarField deltas = edgeNormals() & faPatch::delta(); + const scalarField deltas(edgeNormals() & faPatch::delta()); label sizeby2 = deltas.size()/2; scalar maxMatchError = 0; @@ -213,7 +213,7 @@ void cyclicFaPatch::makeWeights(scalarField& w) const // Make patch edge - neighbour cell distances void cyclicFaPatch::makeDeltaCoeffs(scalarField& dc) const { - scalarField deltas = edgeNormals() & faPatch::delta(); + const scalarField deltas(edgeNormals() & faPatch::delta()); label sizeby2 = deltas.size()/2; for (label edgei = 0; edgei < sizeby2; edgei++) @@ -255,7 +255,7 @@ void Foam::cyclicFaPatch::movePoints(const pointField& p) // Return delta (P to N) vectors across coupled patch tmp cyclicFaPatch::delta() const { - vectorField patchD = faPatch::delta(); + const vectorField patchD(faPatch::delta()); label sizeby2 = patchD.size()/2; tmp tpdv(new vectorField(patchD.size())); diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C index 5061aecb5f..31c7d98a6c 100644 --- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C +++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C @@ -350,7 +350,7 @@ Foam::tmp Foam::faPatch::ngbPolyPatchPointNormals() const tmp tpN(new vectorField(pntEdges.size(), vector::zero)); vectorField& pN = tpN.ref(); - vectorField faceNormals = ngbPolyPatchFaceNormals(); + const vectorField faceNormals(ngbPolyPatchFaceNormals()); forAll(pN, pointI) { diff --git a/src/finiteArea/fields/faPatchFields/basic/basicSymmetry/basicSymmetryFaPatchField.C b/src/finiteArea/fields/faPatchFields/basic/basicSymmetry/basicSymmetryFaPatchField.C index 5a28e331d5..233697a700 100644 --- a/src/finiteArea/fields/faPatchFields/basic/basicSymmetry/basicSymmetryFaPatchField.C +++ b/src/finiteArea/fields/faPatchFields/basic/basicSymmetry/basicSymmetryFaPatchField.C @@ -100,7 +100,8 @@ basicSymmetryFaPatchField::basicSymmetryFaPatchField template tmp > basicSymmetryFaPatchField::snGrad() const { - vectorField nHat = this->patch().edgeNormals(); + const vectorField nHat(this->patch().edgeNormals()); + return ( transform(I - 2.0*sqr(nHat), this->patchInternalField()) @@ -118,7 +119,7 @@ void basicSymmetryFaPatchField::evaluate(const Pstream::commsTypes) this->updateCoeffs(); } - vectorField nHat = this->patch().edgeNormals(); + const vectorField nHat(this->patch().edgeNormals()); Field::operator= ( ( @@ -135,7 +136,7 @@ void basicSymmetryFaPatchField::evaluate(const Pstream::commsTypes) template tmp > basicSymmetryFaPatchField::snGradTransformDiag() const { - vectorField nHat = this->patch().edgeNormals(); + const vectorField nHat(this->patch().edgeNormals()); vectorField diag(nHat.size()); diag.replace(vector::X, mag(nHat.component(vector::X))); diff --git a/src/finiteArea/fields/faPatchFields/constraint/wedge/wedgeFaPatchField.C b/src/finiteArea/fields/faPatchFields/constraint/wedge/wedgeFaPatchField.C index 27a96d83ef..aad90fdf67 100644 --- a/src/finiteArea/fields/faPatchFields/constraint/wedge/wedgeFaPatchField.C +++ b/src/finiteArea/fields/faPatchFields/constraint/wedge/wedgeFaPatchField.C @@ -127,7 +127,7 @@ wedgeFaPatchField::wedgeFaPatchField template tmp > wedgeFaPatchField::snGrad() const { - Field pif = this->patchInternalField(); + const Field pif(this->patchInternalField()); return ( diff --git a/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.C b/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.C index 09d8a57612..6fb0e569f5 100644 --- a/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.C +++ b/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.C @@ -96,9 +96,11 @@ gaussConvectionScheme::famDiv // } // Non-euclidian and other corrections - GeometricField convFluxCorr = + GeometricField convFluxCorr + ( flux(faceFlux, vf) - - faceFlux*tinterpScheme_().euclidianInterpolate(vf);; + - faceFlux*tinterpScheme_().euclidianInterpolate(vf) + ); fam += fac::edgeIntegrate(convFluxCorr); diff --git a/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C b/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C index 0addaac77e..fa67d044f1 100644 --- a/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C +++ b/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C @@ -197,7 +197,8 @@ boundedBackwardFaDdtScheme::facDdt // Calculate unboundedness indicator // Note: all times moved by one because access to internal field // copies current field into the old-time level. - areaScalarField phict = + areaScalarField phict + ( mag ( vf.oldTime().oldTime() @@ -210,13 +211,14 @@ boundedBackwardFaDdtScheme::facDdt - vf.oldTime().oldTime() ) + dimensionedScalar("small", vf.dimensions(), SMALL) - ); + ) + ); - areaScalarField limiter = pos(phict) - pos(phict - scalar(1)); + areaScalarField limiter(pos(phict) - pos(phict - scalar(1))); - areaScalarField coefft = scalar(1) + limiter*deltaT/(deltaT + deltaT0); - areaScalarField coefft00 = limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)); - areaScalarField coefft0 = coefft + coefft00; + areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0)); + areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))); + areaScalarField coefft0(coefft + coefft00); if (mesh().moving()) { @@ -293,7 +295,8 @@ boundedBackwardFaDdtScheme::facDdt0 // Calculate unboundedness indicator // Note: all times moved by one because access to internal field // copies current field into the old-time level. - areaScalarField phict = + areaScalarField phict + ( mag ( vf.oldTime().oldTime() @@ -306,13 +309,14 @@ boundedBackwardFaDdtScheme::facDdt0 - vf.oldTime().oldTime() ) + dimensionedScalar("small", vf.dimensions(), SMALL) - ); + ) + ); - areaScalarField limiter = pos(phict) - pos(phict - scalar(1)); + areaScalarField limiter(pos(phict) - pos(phict - scalar(1))); - areaScalarField coefft = scalar(1) + limiter*deltaT/(deltaT + deltaT0); - areaScalarField coefft00 = limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)); - areaScalarField coefft0 = coefft + coefft00; + areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0)); + areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))); + areaScalarField coefft0(coefft + coefft00); if (mesh().moving()) { @@ -387,7 +391,8 @@ boundedBackwardFaDdtScheme::facDdt // Calculate unboundedness indicator // Note: all times moved by one because access to internal field // copies current field into the old-time level. - areaScalarField phict = + areaScalarField phict + ( mag ( vf.oldTime().oldTime() @@ -400,13 +405,14 @@ boundedBackwardFaDdtScheme::facDdt - vf.oldTime().oldTime() ) + dimensionedScalar("small", vf.dimensions(), SMALL) - ); + ) + ); - areaScalarField limiter = pos(phict) - pos(phict - scalar(1)); + areaScalarField limiter(pos(phict) - pos(phict - scalar(1))); - areaScalarField coefft = scalar(1) + limiter*deltaT/(deltaT + deltaT0); - areaScalarField coefft00 = limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)); - areaScalarField coefft0 = coefft + coefft00; + areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0)); + areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))); + areaScalarField coefft0(coefft + coefft00); if (mesh().moving()) { @@ -483,7 +489,8 @@ boundedBackwardFaDdtScheme::facDdt0 // Calculate unboundedness indicator // Note: all times moved by one because access to internal field // copies current field into the old-time level. - areaScalarField phict = + areaScalarField phict + ( mag ( vf.oldTime().oldTime() @@ -496,13 +503,14 @@ boundedBackwardFaDdtScheme::facDdt0 - vf.oldTime().oldTime() ) + dimensionedScalar("small", vf.dimensions(), SMALL) - ); + ) + ); - areaScalarField limiter = pos(phict) - pos(phict - scalar(1)); + areaScalarField limiter(pos(phict) - pos(phict - scalar(1))); - areaScalarField coefft = scalar(1) + limiter*deltaT/(deltaT + deltaT0); - areaScalarField coefft00 = limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)); - areaScalarField coefft0 = coefft + coefft00; + areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0)); + areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))); + areaScalarField coefft0(coefft + coefft00); if (mesh().moving()) { @@ -577,7 +585,8 @@ boundedBackwardFaDdtScheme::facDdt // Calculate unboundedness indicator // Note: all times moved by one because access to internal field // copies current field into the old-time level. - areaScalarField phict = + areaScalarField phict + ( mag ( rho.oldTime().oldTime()*vf.oldTime().oldTime() @@ -590,13 +599,14 @@ boundedBackwardFaDdtScheme::facDdt - rho.oldTime().oldTime()*vf.oldTime().oldTime() ) + dimensionedScalar("small", rho.dimensions()*vf.dimensions(), SMALL) - ); + ) + ); - areaScalarField limiter = pos(phict) - pos(phict - scalar(1)); + areaScalarField limiter(pos(phict) - pos(phict - scalar(1))); - areaScalarField coefft = scalar(1) + limiter*deltaT/(deltaT + deltaT0); - areaScalarField coefft00 = limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)); - areaScalarField coefft0 = coefft + coefft00; + areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0)); + areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))); + areaScalarField coefft0(coefft + coefft00); if (mesh().moving()) { @@ -677,7 +687,8 @@ boundedBackwardFaDdtScheme::facDdt0 // Calculate unboundedness indicator // Note: all times moved by one because access to internal field // copies current field into the old-time level. - areaScalarField phict = + areaScalarField phict + ( mag ( rho.oldTime().oldTime()*vf.oldTime().oldTime() @@ -690,13 +701,14 @@ boundedBackwardFaDdtScheme::facDdt0 - rho.oldTime().oldTime()*vf.oldTime().oldTime() ) + dimensionedScalar("small", rho.dimensions()*vf.dimensions(), SMALL) - ); + ) + ); - areaScalarField limiter = pos(phict) - pos(phict - scalar(1)); + areaScalarField limiter(pos(phict) - pos(phict - scalar(1))); - areaScalarField coefft = scalar(1) + limiter*deltaT/(deltaT + deltaT0); - areaScalarField coefft00 = limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)); - areaScalarField coefft0 = coefft + coefft00; + areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0)); + areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))); + areaScalarField coefft0(coefft + coefft00); if (mesh().moving()) { @@ -775,7 +787,8 @@ boundedBackwardFaDdtScheme::famDdt // Calculate unboundedness indicator // Note: all times moved by one because access to internal field // copies current field into the old-time level. - scalarField phict = + scalarField phict + ( mag ( vf.oldTime().oldTime().internalField() @@ -788,13 +801,14 @@ boundedBackwardFaDdtScheme::famDdt - vf.oldTime().oldTime().internalField() ) + SMALL - ); + ) + ); scalarField limiter(pos(phict) - pos(phict - 1.0)); - scalarField coefft = 1.0 + limiter*deltaT/(deltaT + deltaT0); - scalarField coefft00 = limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)); - scalarField coefft0 = coefft + coefft00; + scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0)); + scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))); + scalarField coefft0(coefft + coefft00); fam.diag() = (coefft*rDeltaT)*mesh().S(); @@ -845,7 +859,8 @@ boundedBackwardFaDdtScheme::famDdt // Calculate unboundedness indicator // Note: all times moved by one because access to internal field // copies current field into the old-time level. - scalarField phict = + scalarField phict + ( mag ( vf.oldTime().oldTime().internalField() @@ -858,13 +873,14 @@ boundedBackwardFaDdtScheme::famDdt - vf.oldTime().oldTime().internalField() ) + SMALL - ); + ) + ); scalarField limiter(pos(phict) - pos(phict - 1.0)); - scalarField coefft = 1.0 + limiter*deltaT/(deltaT + deltaT0); - scalarField coefft00 = limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)); - scalarField coefft0 = coefft + coefft00; + scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0)); + scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))); + scalarField coefft0(coefft + coefft00); fam.diag() = (coefft*rDeltaT*rho.value())*mesh().S(); @@ -915,7 +931,8 @@ boundedBackwardFaDdtScheme::famDdt // Calculate unboundedness indicator // Note: all times moved by one because access to internal field // copies current field into the old-time level. - scalarField phict = + scalarField phict + ( mag ( rho.oldTime().oldTime().internalField()* @@ -932,13 +949,14 @@ boundedBackwardFaDdtScheme::famDdt vf.oldTime().oldTime().internalField() ) + SMALL - ); + ) + ); scalarField limiter(pos(phict) - pos(phict - 1.0)); - scalarField coefft = 1.0 + limiter*deltaT/(deltaT + deltaT0); - scalarField coefft00 = limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)); - scalarField coefft0 = coefft + coefft00; + scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0)); + scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))); + scalarField coefft0(coefft + coefft00); fam.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().S(); diff --git a/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.C b/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.C index ed8f56d5b4..690634e709 100644 --- a/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.C +++ b/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.C @@ -92,9 +92,11 @@ void gaussGrad::correctBoundaryConditions { if (!vsf.boundaryField()[patchI].coupled()) { - vectorField m = + const vectorField m + ( vsf.mesh().Le().boundaryField()[patchI]/ - vsf.mesh().magLe().boundaryField()[patchI]; + vsf.mesh().magLe().boundaryField()[patchI] + ); gGrad.boundaryFieldRef()[patchI] += m* ( diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C index 7bac7a8e1b..6ed050c518 100644 --- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C +++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C @@ -98,8 +98,8 @@ leastSquaresFaGrad::grad forAll(own, edgei) { - register label ownEdgeI = own[edgei]; - register label neiEdgeI = nei[edgei]; + label ownEdgeI = own[edgei]; + label neiEdgeI = nei[edgei]; Type deltaVsf = vsf[neiEdgeI] - vsf[ownEdgeI]; @@ -117,8 +117,10 @@ leastSquaresFaGrad::grad if (vsf.boundaryField()[patchi].coupled()) { - Field neiVsf = - vsf.boundaryField()[patchi].patchNeighbourField(); + Field neiVsf + ( + vsf.boundaryField()[patchi].patchNeighbourField() + ); forAll(neiVsf, patchEdgeI) { diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C index c914cf187e..c1051741bc 100644 --- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C +++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C @@ -42,7 +42,7 @@ namespace Foam Foam::leastSquaresFaVectors::leastSquaresFaVectors(const faMesh& mesh) : - MeshObject(mesh), + MeshObject(mesh), pVectorsPtr_(NULL), nVectorsPtr_(NULL) {} @@ -135,7 +135,7 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const // HJ, reconsider deltas at the boundary, consistent with FVM // Current implementation is good for fixedValue boudaries, but may // cause problems with fixedGradient. HJ, 4/Oct/2010 - vectorField pd = p.delta(); + const vectorField pd(p.delta()); if (p.coupled()) { @@ -161,7 +161,7 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const // Invert the dd tensor - symmTensorField invDd = inv(dd); + const symmTensorField invDd(inv(dd)); // Revisit all faces and calculate the lsP and lsN vectors @@ -187,7 +187,7 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const const labelUList& edgeFaces = p.edgeFaces(); // Build the d-vectors - vectorField pd = p.delta(); + const vectorField pd(p.delta()); if (p.coupled()) { diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.H b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.H index 1395c8e6eb..fdbd8f5760 100644 --- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.H +++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.H @@ -58,7 +58,7 @@ class mapPolyMesh; class leastSquaresFaVectors : - public MeshObject + public MeshObject { // Private data diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C index faa83de8de..df98e0c9ec 100644 --- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C +++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C @@ -141,7 +141,7 @@ tmp edgeLimitedGrad::grad if (psf.coupled()) { - scalarField psfNei = psf.patchNeighbourField(); + const scalarField psfNei(psf.patchNeighbourField()); forAll(pOwner, pEdgei) { @@ -290,7 +290,7 @@ tmp edgeLimitedGrad::grad if (psf.coupled()) { - vectorField psfNei = psf.patchNeighbourField(); + const vectorField psfNei(psf.patchNeighbourField()); forAll(pOwner, pEdgei) { diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C index 9d52bd87b2..a998b7fa4b 100644 --- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C +++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C @@ -78,7 +78,7 @@ inline void faceLimitedGrad::limitEdge const Type& extrapolate ) { - for(direction cmpt = 0; cmpt < Type::nComponents; cmpt++) + for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++) { faceLimitedGrad::limitEdge ( @@ -145,7 +145,7 @@ tmp faceLimitedGrad::grad if (psf.coupled()) { - scalarField psfNei = psf.patchNeighbourField(); + scalarField psfNei(psf.patchNeighbourField()); forAll(pOwner, pFacei) { @@ -174,7 +174,7 @@ tmp faceLimitedGrad::grad if (k_ < 1.0) { - scalarField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf); + scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf)); maxVsf += maxMinVsf; minVsf -= maxMinVsf; } @@ -293,7 +293,7 @@ tmp faceLimitedGrad::grad if (psf.coupled()) { - vectorField psfNei = psf.patchNeighbourField(); + vectorField psfNei(psf.patchNeighbourField()); forAll(pOwner, pFacei) { @@ -322,7 +322,7 @@ tmp faceLimitedGrad::grad if (k_ < 1.0) { - vectorField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf); + vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf)); maxVsf += maxMinVsf; minVsf -= maxMinVsf; diff --git a/src/finiteArea/finiteArea/laplacianSchemes/gaussFaLaplacianScheme/gaussFaLaplacianScheme.C b/src/finiteArea/finiteArea/laplacianSchemes/gaussFaLaplacianScheme/gaussFaLaplacianScheme.C index 4ed91cf61f..8d1721c1cc 100644 --- a/src/finiteArea/finiteArea/laplacianSchemes/gaussFaLaplacianScheme/gaussFaLaplacianScheme.C +++ b/src/finiteArea/finiteArea/laplacianSchemes/gaussFaLaplacianScheme/gaussFaLaplacianScheme.C @@ -52,7 +52,7 @@ gaussLaplacianScheme::famLaplacian tmp tdeltaCoeffs = this->tlnGradScheme_().deltaCoeffs(vf); const edgeScalarField& deltaCoeffs = tdeltaCoeffs(); - edgeScalarField gammaMagSf = gamma*this->mesh().magLe(); + const edgeScalarField gammaMagSf(gamma*this->mesh().magLe()); tmp > tfam ( diff --git a/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.C b/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.C index 7b2c1c2588..fa3d976e52 100644 --- a/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.C +++ b/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.C @@ -79,7 +79,7 @@ fourthLnGrad::correction ); GeometricField& corr = tcorr.ref(); - edgeVectorField m = mesh.Le()/mesh.magLe(); + edgeVectorField m(mesh.Le()/mesh.magLe()); for (direction cmpt = 0; cmpt < pTraits::nComponents; cmpt++) { diff --git a/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.C b/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.C index 95fd607eb5..a4b969cd9b 100644 --- a/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.C +++ b/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.C @@ -67,10 +67,12 @@ limitedLnGrad::correction const GeometricField& vf ) const { - GeometricField corr = - correctedLnGrad(this->mesh()).correction(vf); + const GeometricField corr + ( + correctedLnGrad(this->mesh()).correction(vf) + ); - edgeScalarField limiter + const edgeScalarField limiter ( min ( diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C index f39499b452..18bee38db4 100644 --- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C +++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C @@ -569,7 +569,7 @@ void edgeInterpolation::makeCorrectionVectors() const if (owner.size() > 0) { - scalarField sinAlpha = deltaCoeffs*mag(CorrVecs.internalField()); + scalarField sinAlpha(deltaCoeffs*mag(CorrVecs.internalField())); forAll(sinAlpha, edgeI) { @@ -660,10 +660,13 @@ void edgeInterpolation::makeSkewCorrectionVectors() const } + + edgeVectorField::Boundary& bSkewCorrVecs = + SkewCorrVecs.boundaryFieldRef(); + forAll(SkewCorrVecs.boundaryField(), patchI) { - faePatchVectorField& patchSkewCorrVecs = - SkewCorrVecs.boundaryFieldRef()[patchI]; + faePatchVectorField& patchSkewCorrVecs = bSkewCorrVecs[patchI]; if (patchSkewCorrVecs.coupled()) { @@ -673,8 +676,7 @@ void edgeInterpolation::makeSkewCorrectionVectors() const const edgeList::subList patchEdges = mesh().boundary()[patchI].patchSlice(edges); - vectorField ngbC = - C.boundaryField()[patchI].patchNeighbourField(); + vectorField ngbC(C.boundaryField()[patchI].patchNeighbourField()); forAll (patchSkewCorrVecs, edgeI) { @@ -683,8 +685,9 @@ void edgeInterpolation::makeSkewCorrectionVectors() const vector S = points[patchEdges[edgeI].start()]; vector e = patchEdges[edgeI].vec(points); - scalar alpha = - ( ( (N - P)^(S - P) )&( (N - P)^e ) )/ - ( ( (N - P)^e )&( (N - P)^e ) ); + scalar alpha = + - (((N - P)^(S - P))&((N - P)^e)) + /(((N - P)^e)&((N - P)^e)); vector E = S + alpha*e; diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C index dbbdd6c1b1..21e1cb8ee1 100644 --- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C +++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C @@ -209,7 +209,7 @@ edgeInterpolationScheme::interpolate Field& sfi = sf.primitiveFieldRef(); - for (register label fi=0; fi::interpolate Field& sfi = sf.primitiveFieldRef(); - for (register label eI = 0; eI < P.size(); eI++) + for (label eI = 0; eI < P.size(); eI++) { // ZT, 22/Apr/2003 const tensorField& curT = mesh.edgeTransformTensors()[eI]; @@ -360,8 +360,8 @@ edgeInterpolationScheme::interpolate label size = vf.boundaryField()[pi].patch().size(); label start = vf.boundaryField()[pi].patch().start(); - Field pOwnVf = vf.boundaryField()[pi].patchInternalField(); - Field pNgbVf = vf.boundaryField()[pi].patchNeighbourField(); + Field pOwnVf(vf.boundaryField()[pi].patchInternalField()); + Field pNgbVf(vf.boundaryField()[pi].patchNeighbourField()); Field& pSf = sf.boundaryFieldRef()[pi]; @@ -449,7 +449,7 @@ edgeInterpolationScheme::euclidianInterpolate Field& sfi = sf.primitiveFieldRef(); - for (register label eI = 0; eI < P.size(); eI++) + for (label eI = 0; eI < P.size(); eI++) { sfi[eI] = lambda[eI]*vfi[P[eI]] + (1 - lambda[eI])*vfi[N[eI]]; } diff --git a/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C b/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C index bf0c11d37b..0a211fcdb2 100644 --- a/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C +++ b/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C @@ -55,11 +55,13 @@ inline tmp limiter(const areaScalarField& phi) return phi; } + inline tmp limiter(const areaVectorField& phi) { return magSqr(phi); } + inline tmp limiter(const areaTensorField& phi) { return magSqr(phi); @@ -88,11 +90,13 @@ tmp faNVDscheme::weights tmp tvf = limiter(phi); const areaScalarField& vf = tvf(); - areaVectorField gradc(fac::grad(vf)); + const areaVectorField gradc(fac::grad(vf)); -// edgeVectorField d = +// edgeVectorField d +// ( // mesh.Le() -// /(mesh.magLe()*mesh.edgeInterpolation::deltaCoeffs()); +// /(mesh.magLe()*mesh.edgeInterpolation::deltaCoeffs()) +// ); // if (!mesh.orthogonal()) // { @@ -110,7 +114,7 @@ tmp faNVDscheme::weights { vector d = vector::zero; - if(edgeFlux_[edge] > 0) + if (edgeFlux_[edge] > 0) { d = c[neighbour[edge]] - c[owner[edge]]; d -= n[owner[edge]]*(n[owner[edge]]&d); @@ -148,41 +152,53 @@ tmp faNVDscheme::weights const scalarField& pEdgeFlux = edgeFlux_.boundaryField()[patchI]; - scalarField pVfP = - vf.boundaryField()[patchI].patchInternalField(); + scalarField pVfP(vf.boundaryField()[patchI].patchInternalField()); - scalarField pVfN = - vf.boundaryField()[patchI].patchNeighbourField(); + scalarField pVfN(vf.boundaryField()[patchI].patchNeighbourField()); - vectorField pGradcP = - gradc.boundaryField()[patchI].patchInternalField(); + vectorField pGradcP + ( + gradc.boundaryField()[patchI].patchInternalField() + ); - vectorField pGradcN = - gradc.boundaryField()[patchI].patchNeighbourField(); + vectorField pGradcN + ( + gradc.boundaryField()[patchI].patchNeighbourField() + ); - vectorField CP = mesh.areaCentres().boundaryField()[patchI] - .patchInternalField(); + vectorField CP + ( + mesh.areaCentres().boundaryField()[patchI].patchInternalField() + ); - vectorField CN = + vectorField CN + ( mesh.areaCentres().boundaryField()[patchI] - .patchNeighbourField(); + .patchNeighbourField() + ); - vectorField nP = + vectorField nP + ( mesh.faceAreaNormals().boundaryField()[patchI] - .patchInternalField(); + .patchInternalField() + ); - vectorField nN = + vectorField nN + ( mesh.faceAreaNormals().boundaryField()[patchI] - .patchNeighbourField(); + .patchNeighbourField() + ); - scalarField pLPN = - mesh.edgeInterpolation::lPN().boundaryField()[patchI]; + scalarField pLPN + ( + mesh.edgeInterpolation::lPN().boundaryField()[patchI] + ); forAll(pWeights, edgeI) { vector d = vector::zero; - if(pEdgeFlux[edgeI] > 0) + if (pEdgeFlux[edgeI] > 0) { d = CN[edgeI] - CP[edgeI]; d -= nP[edgeI]*(nP[edgeI]&d);