diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C index 8eabdce99f..0733802b08 100644 --- a/applications/solvers/DNS/dnsFoam/dnsFoam.C +++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C @@ -86,23 +86,28 @@ int main(int argc, char *argv[]) for (int corr=1; corr<=1; corr++) { volScalarField rAU(1.0/UEqn.A()); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); - U = rAU*UEqn.H(); - phi = (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, U, phi); + surfaceScalarField phiHbyA + ( + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, U, phi) + ); fvScalarMatrix pEqn ( - fvm::laplacian(rAU, p) == fvc::div(phi) + fvm::laplacian(rAU, p) == fvc::div(phiHbyA) ); pEqn.solve(); - phi -= pEqn.flux(); + phi = phiHbyA - pEqn.flux(); #include "continuityErrs.H" - U -= rAU*fvc::grad(p); + U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); } diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H index cb25e7029f..723df7af98 100644 --- a/applications/solvers/combustion/XiFoam/pEqn.H +++ b/applications/solvers/combustion/XiFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -34,19 +35,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -54,7 +58,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -62,7 +66,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H index d1b7135649..00603dc7d9 100644 --- a/applications/solvers/combustion/engineFoam/pEqn.H +++ b/applications/solvers/combustion/engineFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -9,7 +10,7 @@ if (pimple.transonic()) ( "phid", fvc::interpolate(psi) - *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)) + *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U)) ); while (pimple.correctNonOrthogonal()) @@ -31,15 +32,19 @@ if (pimple.transonic()) } else { - phi = fvc::interpolate(rho) - *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho) + *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U)) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -47,7 +52,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -55,7 +60,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H index 872a0513a3..f5364ae314 100644 --- a/applications/solvers/combustion/fireFoam/pEqn.H +++ b/applications/solvers/combustion/fireFoam/pEqn.H @@ -2,25 +2,29 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU)); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); -surfaceScalarField phiU +surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + +surfaceScalarField phiHbyA ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) + - phig ); -phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf(); while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( fvc::ddt(psi, rho)*gh - + fvc::div(phi) + + fvc::div(phiHbyA) + fvm::ddt(psi, p_rgh) - fvm::laplacian(rhorAUf, p_rgh) == @@ -32,7 +36,9 @@ while (pimple.correctNonOrthogonal()) if (pimple.finalNonOrthogonalIter()) { - phi += p_rghEqn.flux(); + phi = phiHbyA + p_rghEqn.flux(); + U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf); + U.correctBoundaryConditions(); } } @@ -41,8 +47,6 @@ p = p_rgh + rho*gh; #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf); -U.correctBoundaryConditions(); K = 0.5*magSqr(U); dpdt = fvc::ddt(p); diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H index cb25e7029f..723df7af98 100644 --- a/applications/solvers/combustion/reactingFoam/pEqn.H +++ b/applications/solvers/combustion/reactingFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -34,19 +35,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -54,7 +58,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -62,7 +66,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/combustion/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/rhoReactingFoam/pEqn.H index 4e94de243c..ef11c677ad 100644 --- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H +++ b/applications/solvers/combustion/rhoReactingFoam/pEqn.H @@ -6,27 +6,25 @@ thermo.rho() -= psi*p; volScalarField rAU(1.0/UEqn.A()); - U = rAU*UEqn.H(); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); if (pimple.transonic()) { - surfaceScalarField phiv + surfaceScalarField phiHbyA ( - (fvc::interpolate(U) & mesh.Sf()) + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - phi = fvc::interpolate(rho)*phiv; + surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA); - surfaceScalarField phid - ( - "phid", - fvc::interpolate(thermo.psi())*phiv - ); + phiHbyA *= fvc::interpolate(rho); fvScalarMatrix pDDtEqn ( - fvc::ddt(rho) + fvc::div(phi) + fvc::ddt(rho) + fvc::div(phiHbyA) + correction(fvm::ddt(psi, p) + fvm::div(phid, p)) ); @@ -42,23 +40,26 @@ if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); fvScalarMatrix pDDtEqn ( fvc::ddt(rho) + psi*correction(fvm::ddt(p)) - + fvc::div(phi) + + fvc::div(phiHbyA) ); while (pimple.correctNonOrthogonal()) @@ -73,7 +74,7 @@ if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -84,7 +85,7 @@ #include "rhoEqn.H" #include "compressibleContinuityErrs.H" - U -= rAU*fvc::grad(p); + U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H index f7bb932fdc..cfe34f2873 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H @@ -3,7 +3,8 @@ rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); -U = rAU*UEqn().H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn().H(); if (pimple.nCorrPISO() <= 1) { @@ -17,7 +18,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -41,12 +42,15 @@ if (pimple.transonic()) } else { - phi = - fvc::interpolate(rho)* - ( - (fvc::interpolate(U) & mesh.Sf()) + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { @@ -54,7 +58,7 @@ else fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -62,7 +66,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -81,7 +85,7 @@ rho.relax(); Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H index 831f139e81..c5c7602a43 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H @@ -4,7 +4,8 @@ rho = min(rho, rhoMax); rho.relax(); volScalarField rAU(1.0/UEqn().A()); -U = rAU*UEqn().H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn().H(); if (pimple.nCorrPISO() <= 1) { @@ -18,7 +19,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -43,13 +44,17 @@ if (pimple.transonic()) } else { - phi = - fvc::interpolate(rho)* - ( - (fvc::interpolate(U) & mesh.Sf()) + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); - mrfZones.relativeFlux(fvc::interpolate(rho), phi); + ) + ); + + mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA); while (pimple.correctNonOrthogonal()) { @@ -57,7 +62,7 @@ else fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) ); @@ -65,7 +70,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -83,7 +88,7 @@ rho.relax(); Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index 4d3f1ac76e..56c444cdab 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -4,7 +4,9 @@ rho = min(rho, rhoMax); rho.relax(); volScalarField rAU(1.0/UEqn().A()); -U = rAU*UEqn().H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn().H(); + UEqn.clear(); bool closedVolume = false; @@ -14,7 +16,7 @@ if (simple.transonic()) surfaceScalarField phid ( "phid", - fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf()) + fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf()) ); while (simple.correctNonOrthogonal()) @@ -40,14 +42,19 @@ if (simple.transonic()) } else { - phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); - closedVolume = adjustPhi(phi, U, p); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); + + closedVolume = adjustPhi(phiHbyA, U, p); while (simple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( - fvm::laplacian(rho*rAU, p) == fvc::div(phi) + fvm::laplacian(rho*rAU, p) == fvc::div(phiHbyA) ); pEqn.setReference(pRefCell, pRefValue); @@ -56,7 +63,7 @@ else if (simple.finalNonOrthogonalIter()) { - phi -= pEqn.flux(); + phi = phiHbyA - pEqn.flux(); } } } @@ -67,7 +74,7 @@ else // Explicitly relax pressure for momentum corrector p.relax(); -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); // For closed-volume cases adjust the pressure and density levels diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H index b1f83304cd..c3e0a43a15 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H @@ -1,10 +1,12 @@ +volVectorField HbyA("HbyA", U); + if (pressureImplicitPorosity) { - U = trTU()&UEqn().H(); + HbyA = trTU() & UEqn().H(); } else { - U = trAU()*UEqn().H(); + HbyA = trAU()*UEqn().H(); } UEqn.clear(); @@ -16,7 +18,7 @@ if (simple.transonic()) surfaceScalarField phid ( "phid", - fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf()) + fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf()) ); mrfZones.relativeFlux(fvc::interpolate(psi), phid); @@ -45,10 +47,15 @@ if (simple.transonic()) } else { - phi = fvc::interpolate(rho*U) & mesh.Sf(); - mrfZones.relativeFlux(fvc::interpolate(rho), phi); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho*HbyA) & mesh.Sf() + ); - closedVolume = adjustPhi(phi, U, p); + mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA); + + closedVolume = adjustPhi(phiHbyA, U, p); while (simple.correctNonOrthogonal()) { @@ -56,11 +63,11 @@ else if (pressureImplicitPorosity) { - tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phi)); + tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phiHbyA)); } else { - tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phi)); + tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phiHbyA)); } tpEqn().setReference(pRefCell, pRefValue); @@ -69,7 +76,7 @@ else if (simple.finalNonOrthogonalIter()) { - phi -= tpEqn().flux(); + phi = phiHbyA - tpEqn().flux(); } } } @@ -81,11 +88,11 @@ p.relax(); if (pressureImplicitPorosity) { - U -= trTU()&fvc::grad(p); + U = HbyA - (trTU() & fvc::grad(p)); } else { - U -= trAU()*fvc::grad(p); + U = HbyA - trAU()*fvc::grad(p); } U.correctBoundaryConditions(); diff --git a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H index ccc7a1b21e..a4d9325522 100644 --- a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H @@ -7,7 +7,10 @@ volScalarField p0(p); volScalarField AU(UEqn().A()); volScalarField AtU(AU - UEqn().H1()); -U = UEqn().H()/AU; + +volVectorField HbyA("HbyA", U); +HbyA = UEqn().H()/AU; + UEqn.clear(); bool closedVolume = false; @@ -19,7 +22,7 @@ if (simple.transonic()) surfaceScalarField phid ( "phid", - fvc::interpolate(psi*U) & mesh.Sf() + fvc::interpolate(psi*HbyA) & mesh.Sf() ); surfaceScalarField phic @@ -56,13 +59,18 @@ else { while (simple.correctNonOrthogonal()) { - phi = fvc::interpolate(rho*U) & mesh.Sf(); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho*HbyA) & mesh.Sf() + ); + closedVolume = adjustPhi(phi, U, p); phi += fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf(); fvScalarMatrix pEqn ( - fvc::div(phi) + fvc::div(phiHbyA) //- fvm::laplacian(rho/AU, p) - fvm::laplacian(rho/AtU, p) ); @@ -73,7 +81,7 @@ else if (simple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -85,8 +93,8 @@ else // Explicitly relax pressure for momentum corrector p.relax(); -U -= (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU); -//U -= fvc::grad(p)/AU; +U = HbyA - (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU); +//U = HbyA - fvc::grad(p)/AU; U.correctBoundaryConditions(); diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H index a4311dfd2e..46cc36216c 100644 --- a/applications/solvers/compressible/sonicFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/pEqn.H @@ -1,14 +1,15 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); surfaceScalarField phid ( "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -33,5 +34,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H index 98cc22d15a..cad5bb7217 100644 --- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H @@ -1,14 +1,15 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = UEqn.H()/UEqn.A(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); surfaceScalarField phid ( "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U) ) ); @@ -29,5 +30,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C index a91ee8885e..51f04820eb 100644 --- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C +++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C @@ -93,17 +93,21 @@ int main(int argc, char *argv[]) for (int corr=0; corr 0) + surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf()); + if (pZones.size() == 0) { - // ddtPhiCorr not well defined for cases with porosity - phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); - } - else - { - phi = - fvc::interpolate(rho) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + // ddtPhiCorr only used without porosity + phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi); } + phiHbyA *= fvc::interpolate(rho); + + fvScalarMatrix pDDtEqn ( fvc::ddt(rho) + psi*correction(fvm::ddt(p)) - + fvc::div(phi) + + fvc::div(phiHbyA) == parcels.Srho() + sources(psi, p, rho.name()) @@ -46,7 +42,7 @@ if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } @@ -58,7 +54,7 @@ #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent #include "compressibleContinuityErrs.H" - U -= rAU*fvc::grad(p); + U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); sources.correct(U); diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H index cb8553bbfb..3a61b28e60 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*(UEqn == sources(rho, U))().H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*(UEqn == sources(rho, U))().H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -39,19 +40,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) == coalParcels.Srho() @@ -64,7 +68,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -72,7 +76,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); sources.correct(U); diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H index 6ed7affe77..4c9b41bba4 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H @@ -6,27 +6,23 @@ thermo.rho() -= psi*p; volScalarField rAU(1.0/UEqn.A()); - U = rAU*(UEqn == sources(rho, U))().H(); + volVectorField HbyA("HbyA", U); + HbyA = rAU*(UEqn == sources(rho, U))().H(); - if (pZones.size() > 0) + surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf()); + if (pZones.size() == 0) { - // ddtPhiCorr not well defined for cases with porosity - phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); - } - else - { - phi = - fvc::interpolate(rho) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + // ddtPhiCorr only used without porosity + phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi); } + phiHbyA *= fvc::interpolate(rho); + + fvScalarMatrix pDDtEqn ( fvc::ddt(rho) + psi*correction(fvm::ddt(p)) - + fvc::div(phi) + + fvc::div(phiHbyA) == parcels.Srho() + sources(psi, p, rho.name()) @@ -47,7 +43,7 @@ if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } @@ -57,7 +53,7 @@ #include "rhoEqn.H" #include "compressibleContinuityErrs.H" - U -= rAU*fvc::grad(p); + U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); sources.correct(U); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H index 872a0513a3..f5364ae314 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H @@ -2,25 +2,29 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU)); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); -surfaceScalarField phiU +surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + +surfaceScalarField phiHbyA ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) + - phig ); -phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf(); while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( fvc::ddt(psi, rho)*gh - + fvc::div(phi) + + fvc::div(phiHbyA) + fvm::ddt(psi, p_rgh) - fvm::laplacian(rhorAUf, p_rgh) == @@ -32,7 +36,9 @@ while (pimple.correctNonOrthogonal()) if (pimple.finalNonOrthogonalIter()) { - phi += p_rghEqn.flux(); + phi = phiHbyA + p_rghEqn.flux(); + U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf); + U.correctBoundaryConditions(); } } @@ -41,8 +47,6 @@ p = p_rgh + rho*gh; #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf); -U.correctBoundaryConditions(); K = 0.5*magSqr(U); dpdt = fvc::ddt(p); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H index 5ee551a0ea..c504fdf2ea 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -36,19 +37,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - (fvc::interpolate(U) & mesh.Sf()) + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) == parcels.Srho() @@ -58,7 +62,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -66,7 +70,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H index 5109788ea1..5c4fa0adc1 100644 --- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H +++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H @@ -1,7 +1,8 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -U = rAU*UEqn.H(); +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); if (pimple.transonic()) { @@ -10,7 +11,7 @@ if (pimple.transonic()) "phid", fvc::interpolate(psi) *( - ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)) + ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U)) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); @@ -36,19 +37,22 @@ if (pimple.transonic()) } else { - phi = + surfaceScalarField phiHbyA + ( + "phiHbyA", fvc::interpolate(rho) *( - ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)) + ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U)) + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); + ) + ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) - + fvc::div(phi) + + fvc::div(phiHbyA) - fvm::laplacian(rho*rAU, p) == parcels.Srho() @@ -58,7 +62,7 @@ else if (pimple.finalNonOrthogonalIter()) { - phi += pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -66,7 +70,7 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H index 6fa1965a24..34898172ce 100644 --- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H +++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H @@ -9,11 +9,13 @@ )/psi; } - surfaceScalarField rhof(fvc::interpolate(rho, "rhof")); + surfaceScalarField rhof("rhof", fvc::interpolate(rho)); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU)); - volVectorField HbyA("HbyA", rAU*UEqn.H()); + + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); phiv = (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phiv); @@ -22,8 +24,6 @@ phiv -= phiGradp/rhof; - #include "resetPhivPatches.H" - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn diff --git a/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H b/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H deleted file mode 100644 index e7d0c2f93e..0000000000 --- a/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H +++ /dev/null @@ -1,15 +0,0 @@ -fvsPatchScalarFieldField& phiPatches = phi.boundaryField(); -const fvPatchScalarFieldField& rhoPatches = rho.boundaryField(); -const fvPatchVectorFieldField& Upatches = U.boundaryField(); -const fvsPatchVectorFieldField& SfPatches = mesh.Sf().boundaryField(); - -forAll(phiPatches, patchI) -{ - if (phi.boundaryField().types()[patchI] == "calculated") - { - calculatedFvsPatchScalarField& phiPatch = - refCast(phiPatches[patchI]); - - phiPatch == ((rhoPatches[patchI]*Upatches[patchI]) & SfPatches[patchI]); - } -} diff --git a/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H b/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H deleted file mode 100644 index 2f62dcfc37..0000000000 --- a/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H +++ /dev/null @@ -1,19 +0,0 @@ -surfaceScalarField::GeometricBoundaryField& phivPatches = - phiv.boundaryField(); - -const volVectorField::GeometricBoundaryField& Upatches = - U.boundaryField(); - -const surfaceVectorField::GeometricBoundaryField& SfPatches = - mesh.Sf().boundaryField(); - -forAll(phivPatches, patchI) -{ - if (phiv.boundaryField().types()[patchI] == "calculated") - { - calculatedFvsPatchScalarField& phivPatch = - refCast(phivPatches[patchI]); - - phivPatch == (Upatches[patchI] & SfPatches[patchI]); - } -} diff --git a/applications/test/primitivePatch/Test-PrimitivePatch.C b/applications/test/primitivePatch/Test-PrimitivePatch.C index 634b47130b..6f6d998c7d 100644 --- a/applications/test/primitivePatch/Test-PrimitivePatch.C +++ b/applications/test/primitivePatch/Test-PrimitivePatch.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -26,6 +26,7 @@ Description \*---------------------------------------------------------------------------*/ +#include "PrimitivePatch.H" #include "argList.H" #include "Time.H" #include "polyMesh.H" @@ -223,25 +224,40 @@ int main(int argc, char *argv[]) Info<< "Patch:" << patch.name() << endl; - PrimitivePatch pp(patch, patch.points()); + // Test addressing + { + myPrimitivePatch pp(patch, patch.points()); - const pointField& localPoints = pp.localPoints(); - const faceList& localFaces = pp.localFaces(); - const labelListList& faceFaces = pp.faceFaces(); - const edgeList& edges = pp.edges(); - const labelListList& edgeFaces = pp.edgeFaces(); - const labelListList& faceEdges = pp.faceEdges(); + const pointField& localPoints = pp.localPoints(); + const faceList& localFaces = pp.localFaces(); + const labelListList& faceFaces = pp.faceFaces(); + const edgeList& edges = pp.edges(); + const labelListList& edgeFaces = pp.edgeFaces(); + const labelListList& faceEdges = pp.faceEdges(); - checkFaceEdges(localFaces, edges, faceEdges); + checkFaceEdges(localFaces, edges, faceEdges); - writeEdges(localPoints, edges, pp.nInternalEdges()); + writeEdges(localPoints, edges, pp.nInternalEdges()); - writeFaceEdges(localPoints, edges, faceEdges); + writeFaceEdges(localPoints, edges, faceEdges); - writeEdgeFaces(localPoints, localFaces, edges, edgeFaces); + writeEdgeFaces(localPoints, localFaces, edges, edgeFaces); - writeFaceFaces(localPoints, localFaces, faceFaces); + writeFaceFaces(localPoints, localFaces, faceFaces); + } + + // Test construction from Xfer + { + faceList patchFaces = patch; + pointField allPoints = patch.points(); + + PrimitivePatch storedPatch + ( + patchFaces.xfer(), + allPoints.xfer() + ); + } Info<< "End\n" << endl; diff --git a/applications/utilities/mesh/conversion/vtkUnstructuredToFoam/Make/files b/applications/utilities/mesh/conversion/vtkUnstructuredToFoam/Make/files new file mode 100644 index 0000000000..d0a361873c --- /dev/null +++ b/applications/utilities/mesh/conversion/vtkUnstructuredToFoam/Make/files @@ -0,0 +1,4 @@ + +vtkUnstructuredToFoam.C + +EXE = $(FOAM_APPBIN)/vtkUnstructuredToFoam diff --git a/applications/utilities/mesh/conversion/vtkUnstructuredToFoam/Make/options b/applications/utilities/mesh/conversion/vtkUnstructuredToFoam/Make/options new file mode 100644 index 0000000000..7ce182425d --- /dev/null +++ b/applications/utilities/mesh/conversion/vtkUnstructuredToFoam/Make/options @@ -0,0 +1,5 @@ +EXE_INC = \ + -I$(LIB_SRC)/fileFormats/lnInclude + +EXE_LIBS = \ + -lfileFormats diff --git a/applications/utilities/mesh/conversion/vtkUnstructuredToFoam/vtkUnstructuredToFoam.C b/applications/utilities/mesh/conversion/vtkUnstructuredToFoam/vtkUnstructuredToFoam.C new file mode 100644 index 0000000000..753246e43e --- /dev/null +++ b/applications/utilities/mesh/conversion/vtkUnstructuredToFoam/vtkUnstructuredToFoam.C @@ -0,0 +1,87 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Description + Converts ascii .vtk (legacy format) file generated by vtk/paraview. + + Note: the .vtk format does not contain any boundary information. It is + purely a description of the internal mesh. + + Not extensively tested. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Time.H" +#include "polyMesh.H" +#include "IFstream.H" +#include "vtkUnstructuredReader.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Main program: + +int main(int argc, char *argv[]) +{ + argList::noParallel(); + argList::validArgs.append(".vtk ascii file"); + +# include "setRootCase.H" +# include "createTime.H" + + IFstream mshStream(args[1]); + + vtkUnstructuredReader reader(runTime, mshStream); + + polyMesh mesh + ( + IOobject + ( + polyMesh::defaultRegion, + runTime.constant(), + runTime + ), + xferMove(reader.points()), + reader.cells(), + faceListList(0), + wordList(0), + wordList(0), + "defaultFaces", + polyPatch::typeName, + wordList(0) + ); + + Info<< "Writing mesh ..." << endl; + + mesh.write(); + + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C index 7f4c7abb6e..9bc2512f58 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C @@ -41,10 +41,10 @@ void Foam::conformalVoronoiMesh::conformToSurface() { reconformationMode reconfMode = reconformationControl(); - if (Pstream::parRun()) - { - seedProcessorBoundarySurfaces(true); - } +// if (Pstream::parRun()) +// { +// seedProcessorBoundarySurfaces(true); +// } if (reconfMode == rmNone) { @@ -74,15 +74,15 @@ void Foam::conformalVoronoiMesh::conformToSurface() storeSurfaceConformation(); } - if (Pstream::parRun()) - { - label nFarPoints = removeProcessorBoundarySeeds(true); - - reduce(nFarPoints, sumOp