Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/DNS/dnsFoam/dnsFoam.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
@@ -85,7 +85,7 @@ int main(int argc, char *argv[])
 
         for (int corr=1; corr<=1; corr++)
         {
-            volScalarField rAU(1.0/UEqn.A());
+            volScalarField rAU("Dp", 1.0/UEqn.A());
             volVectorField HbyA("HbyA", U);
             HbyA = rAU*UEqn.H();
 

applications/solvers/combustion/chemFoam/createSingleCellMesh.H

@@ -23,7 +23,7 @@ fvMesh mesh
         fvMesh::defaultRegion,
         runTime.timeName(),
         runTime,
-        IOobject::MUST_READ
+        IOobject::READ_IF_PRESENT
     ),
     xferMove<Field<vector> >(points),
     faces.xfer(),

applications/solvers/combustion/fireFoam/pEqn.H

@@ -1,11 +1,11 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
-surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
@@ -15,7 +15,7 @@ surfaceScalarField phiHbyA
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
-  - phig
+  + phig
 );
 
 
@@ -37,7 +37,7 @@ while (pimple.correctNonOrthogonal())
     if (pimple.finalNonOrthogonalIter())
     {
         phi = phiHbyA + p_rghEqn.flux();
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf);
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
         U.correctBoundaryConditions();
     }
 }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -1,18 +1,18 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
+    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, U, phi)
-      - phig
+      + phig
     );
 
     while (pimple.correctNonOrthogonal())
@@ -36,7 +36,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -1,12 +1,12 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn().A());
-    surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();
     UEqn.clear();
 
-    surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
+    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
@@ -16,7 +16,7 @@
 
     adjustPhi(phiHbyA, U, p_rgh);
 
-    phiHbyA -= phig;
+    phiHbyA += phig;
 
     while (simple.correctNonOrthogonal())
     {
@@ -39,7 +39,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -6,7 +6,7 @@
     thermo.rho() -= psi*p_rgh;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -3,13 +3,13 @@
     rho.relax();
 
     volScalarField rAU(1.0/UEqn().A());
-    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();
     UEqn.clear();
 
-    surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
@@ -19,7 +19,7 @@
 
     bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
 
-    phiHbyA -= phig;
+    phiHbyA += phig;
 
     while (simple.correctNonOrthogonal())
     {
@@ -41,7 +41,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -5,7 +5,7 @@
     rho.relax();
 
     volScalarField rAU(1.0/UEqn().A());
-    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 
     U = rAU*UEqn().H();
     UEqn.clear();
@@ -15,8 +15,8 @@
     dimensionedScalar compressibility = fvc::domainIntegrate(psi);
     bool compressible = (compressibility.value() > SMALL);
 
-    surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-    phi -= buoyancyPhi;
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+    phi += phig;
 
     // Solve pressure
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -44,7 +44,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
+            U += rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -6,7 +6,7 @@
     rho = thermo.rho();
 
     volScalarField rAU(1.0/UEqn().A());
-    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();

applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C

@@ -61,6 +61,7 @@ int main(int argc, char *argv[])
             fvm::ddt(U)
           + fvm::div(phi, U)
           - fvm::laplacian(fluid.nu(), U)
+          - (fvc::grad(U) & fvc::grad(fluid.nu()))
         );
 
         solve(UEqn == -fvc::grad(p));

applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H

@@ -1,5 +1,5 @@
 volScalarField rAU(1.0/UEqn().A());
-surfaceScalarField rAUf(rAU.name() + 'f', fvc::interpolate(rAU));
+surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*(UEqn() == sources(U))().H();

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -1,11 +1,11 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
-surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
@@ -15,7 +15,7 @@ surfaceScalarField phiHbyA
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
-  - phig
+  + phig
 );
 
 
@@ -37,7 +37,7 @@ while (pimple.correctNonOrthogonal())
     if (pimple.finalNonOrthogonalIter())
     {
         phi = phiHbyA + p_rghEqn.flux();
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf);
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
         U.correctBoundaryConditions();
     }
 }

applications/solvers/multiphase/bubbleFoam/pEqn.H

@@ -8,16 +8,21 @@
     surfaceScalarField rAU1f(fvc::interpolate(rAU1));
     surfaceScalarField rAU2f(fvc::interpolate(rAU2));
 
-    U1 = rAU1*U1Eqn.H();
-    U2 = rAU2*U2Eqn.H();
+    volVectorField HbyA1("HbyA1", U1);
+    HbyA1 = rAU1*U1Eqn.H();
+
+    volVectorField HbyA2("HbyA2", U2);
+    HbyA2 = rAU2*U2Eqn.H();
 
     surfaceScalarField phiDrag1
     (
-        fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2 + rAU1f*(g & mesh.Sf())
+        fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2
+      + rAU1f*(g & mesh.Sf())
     );
     surfaceScalarField phiDrag2
     (
-        fvc::interpolate(alpha1/rho2*dragCoef*rAU2)*phi1 + rAU2f*(g & mesh.Sf())
+        fvc::interpolate(alpha1/rho2*dragCoef*rAU2)*phi1
+      + rAU2f*(g & mesh.Sf())
     );
 
     forAll(p.boundaryField(), patchi)
@@ -29,16 +34,25 @@
         }
     }
 
-    phi1 = (fvc::interpolate(U1) & mesh.Sf()) + fvc::ddtPhiCorr(rAU1, U1, phi1)
-        + phiDrag1;
-    phi2 = (fvc::interpolate(U2) & mesh.Sf()) + fvc::ddtPhiCorr(rAU2, U2, phi2)
-        + phiDrag2;
+    surfaceScalarField phiHbyA1
+    (
+        (fvc::interpolate(HbyA1) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU1, U1, phi1)
+      + phiDrag1
+    );
 
-    phi = alpha1f*phi1 + alpha2f*phi2;
+    surfaceScalarField phiHbyA2
+    (
+        (fvc::interpolate(HbyA2) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU2, U2, phi2)
+      + phiDrag2
+    );
+
+    surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
     surfaceScalarField Dp
     (
-        "(rho*(1|A(U)))",
+        "Dp",
         alpha1f*rAU1f/rho1 + alpha2f*rAU2f/rho2
     );
 
@@ -46,7 +60,7 @@
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(Dp, p) == fvc::div(phi)
+            fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -57,19 +71,17 @@
         {
             surfaceScalarField SfGradp(pEqn.flux()/Dp);
 
-            phi1 -= rAU1f*SfGradp/rho1;
-            phi2 -= rAU2f*SfGradp/rho2;
+            phi1 = phiHbyA1 - rAU1f*SfGradp/rho1;
+            phi2 = phiHbyA2 - rAU2f*SfGradp/rho2;
             phi = alpha1f*phi1 + alpha2f*phi2;
 
             p.relax();
             SfGradp = pEqn.flux()/Dp;
 
-            U1 += (fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1));
-            //U1 += rAU1*(fvc::reconstruct(phiDrag1/rAU1f - SfGradp/rho1));
+            U1 = HbyA1 + (fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1));
             U1.correctBoundaryConditions();
 
-            U2 += (fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2));
-            //U2 += rAU2*(fvc::reconstruct(phiDrag2/rAU2f - SfGradp/rho2));
+            U2 = HbyA2 + (fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2));
             U2.correctBoundaryConditions();
 
             U = alpha1*U1 + alpha2*U2;

applications/solvers/multiphase/cavitatingFoam/createFields.H

@@ -25,38 +25,6 @@
         mesh
     );
 
-    volScalarField gamma
-    (
-        IOobject
-        (
-            "gamma",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0))
-    );
-    gamma.oldTime();
-
-    Info<< "Creating compressibilityModel\n" << endl;
-    autoPtr<barotropicCompressibilityModel> psiModel =
-        barotropicCompressibilityModel::New
-        (
-            thermodynamicProperties,
-            gamma
-        );
-
-    const volScalarField& psi = psiModel->psi();
-
-    rho == max
-    (
-        psi*p
-      + (1.0 - gamma)*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
-        rhoMin
-    );
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -78,6 +46,27 @@
 
     twoPhaseMixture twoPhaseProperties(U, phiv, "gamma");
 
+    volScalarField& gamma(twoPhaseProperties.alpha1());
+    gamma.oldTime();
+
+    Info<< "Creating compressibilityModel\n" << endl;
+    autoPtr<barotropicCompressibilityModel> psiModel =
+        barotropicCompressibilityModel::New
+        (
+            thermodynamicProperties,
+            gamma
+        );
+
+    const volScalarField& psi = psiModel->psi();
+
+    rho == max
+    (
+        psi*p
+      + (1.0 - gamma)*rhol0
+      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
+        rhoMin
+    );
+
     // Create incompressible turbulence model
     autoPtr<incompressible::turbulenceModel> turbulence
     (

applications/solvers/multiphase/cavitatingFoam/pEqn.H

@@ -12,7 +12,7 @@
     surfaceScalarField rhof("rhof", fvc::interpolate(rho));
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", rhof*fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();

applications/solvers/multiphase/compressibleInterFoam/Allwclean

@@ -2,6 +2,7 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
+wclean libso phaseEquationsOfState
 wclean
 wclean compressibleInterDyMFoam
 

applications/solvers/multiphase/compressibleInterFoam/Allwmake

@@ -2,6 +2,7 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
+wmake libso phaseEquationsOfState
 wmake
 wmake compressibleInterDyMFoam
 

applications/solvers/multiphase/compressibleInterFoam/Make/files

@@ -1,3 +1,4 @@
+derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C
 compressibleInterFoam.C
 
 EXE = $(FOAM_APPBIN)/compressibleInterFoam

applications/solvers/multiphase/compressibleInterFoam/Make/options

@@ -2,12 +2,14 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -IphaseEquationsOfState/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -ltwoPhaseInterfaceProperties \
     -lincompressibleTransportModels \
+    -lphaseEquationsOfState \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleLESModels \

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -0,0 +1,20 @@
+{
+    volScalarField kByCv
+    (
+        "kByCv",
+        (alpha1*k1/Cv1 + alpha2*k2/Cv2)
+      + (alpha1*rho1 + alpha2*rho2)*turbulence->nut()
+    );
+
+    solve
+    (
+        fvm::ddt(rho, T)
+      + fvm::div(rhoPhi, T)
+      - fvm::laplacian(kByCv, T)
+      + p*fvc::div(phi)*(alpha1/Cv1 + alpha2/Cv2)
+    );
+
+    // Update compressibilities
+    psi1 = eos1->psi(p, T);
+    psi2 = eos2->psi(p, T);
+}

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options

@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I../phaseEquationsOfState/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
@@ -12,6 +13,7 @@ EXE_INC = \
 EXE_LIBS = \
     -ltwoPhaseInterfaceProperties \
     -lincompressibleTransportModels \
+    -lphaseEquationsOfState \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleLESModels \

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.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
@@ -25,7 +25,7 @@ Application
     compressibleInterDyMFoam
 
 Description
-    Solver for 2 compressible, isothermal immiscible fluids using a VOF
+    Solver for 2 compressible, non-isothermal immiscible fluids using a VOF
     (volume of fluid) phase-fraction based interface capturing approach,
     with optional mesh motion and mesh topology changes including adaptive
     re-meshing.
@@ -43,6 +43,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
+#include "phaseEquationOfState.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 
@@ -124,11 +125,15 @@ int main(int argc, char *argv[])
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             #include "UEqn.H"
+            #include "TEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())
             {
                 #include "pEqn.H"
+
+                // Make the fluxes relative to the mesh motion
+                fvc::makeRelative(phi, U);
             }
         }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -1,88 +0,0 @@
-{
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
-
-    tmp<fvScalarMatrix> p_rghEqnComp;
-
-    if (pimple.transonic())
-    {
-        p_rghEqnComp =
-        (
-            fvm::ddt(p_rgh)
-          + fvm::div(phi, p_rgh)
-          - fvm::Sp(fvc::div(phi), p_rgh)
-        );
-    }
-    else
-    {
-        p_rghEqnComp =
-        (
-            fvm::ddt(p_rgh)
-          + fvc::div(phi, p_rgh)
-          - fvc::Sp(fvc::div(phi), p_rgh)
-        );
-    }
-
-
-    U = rAU*UEqn.H();
-
-    surfaceScalarField phiU
-    (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
-    );
-
-    phi = phiU +
-        (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix p_rghEqnIncomp
-        (
-            fvc::div(phi)
-          - fvm::laplacian(rAUf, p_rgh)
-        );
-
-        solve
-        (
-            (
-                max(alpha1, scalar(0))*(psi1/rho1)
-              + max(alpha2, scalar(0))*(psi2/rho2)
-            )
-           *p_rghEqnComp()
-          + p_rghEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
-        );
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            dgdt =
-                (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(p_rghEqnComp & p_rgh);
-            phi += p_rghEqnIncomp.flux();
-        }
-    }
-
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
-    p = max
-    (
-        (p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
-       /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-        pMin
-    );
-
-    rho1 = rho10 + psi1*p;
-    rho2 = rho20 + psi2*p;
-
-    Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
-
-    // Make the fluxes relative to the mesh motion
-    fvc::makeRelative(phi, U);
-}

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.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
@@ -25,7 +25,7 @@ Application
     compressibleInterFoam
 
 Description
-    Solver for 2 compressible, isothermal immiscible fluids using a VOF
+    Solver for 2 compressible, non-isothermal immiscible fluids using a VOF
     (volume of fluid) phase-fraction based interface capturing approach.
 
     The momentum and other fluid properties are of the "mixture" and a single
@@ -40,6 +40,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
+#include "phaseEquationOfState.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 
@@ -82,6 +83,7 @@ int main(int argc, char *argv[])
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             #include "UEqn.H"
+            #include "TEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -12,23 +12,6 @@
         mesh
     );
 
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Calculating field alpha1\n" << endl;
-    volScalarField alpha2("alpha2", scalar(1) - alpha1);
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -45,48 +28,20 @@
 
     #include "createPhi.H"
 
-
-    Info<< "Reading transportProperties\n" << endl;
-    twoPhaseMixture twoPhaseProperties(U, phi);
-
-    dimensionedScalar rho10
+    Info<< "Reading field T\n" << endl;
+    volScalarField T
     (
-        twoPhaseProperties.subDict
+        IOobject
         (
-            twoPhaseProperties.phase1Name()
-        ).lookup("rho0")
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
     );
 
-    dimensionedScalar rho20
-    (
-        twoPhaseProperties.subDict
-        (
-            twoPhaseProperties.phase2Name()
-        ).lookup("rho0")
-    );
-
-    dimensionedScalar psi1
-    (
-        twoPhaseProperties.subDict
-        (
-            twoPhaseProperties.phase1Name()
-        ).lookup("psi")
-    );
-
-    dimensionedScalar psi2
-    (
-        twoPhaseProperties.subDict
-        (
-            twoPhaseProperties.phase2Name()
-        ).lookup("psi")
-    );
-
-    dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
-
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
     volScalarField p
     (
         IOobject
@@ -94,19 +49,120 @@
             "p",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        max
+        p_rgh
+    );
+
+
+    Info<< "Reading transportProperties\n" << endl;
+    twoPhaseMixture twoPhaseProperties(U, phi);
+
+    volScalarField& alpha1(twoPhaseProperties.alpha1());
+
+    Info<< "Calculating phase-fraction alpha" << twoPhaseProperties.phase2Name()
+        << nl << endl;
+    volScalarField alpha2
     (
-            (p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
+        "alpha" + twoPhaseProperties.phase2Name(),
+        scalar(1) - alpha1
+    );
+
+    dimensionedScalar k1
+    (
+        "k",
+        dimensionSet(1, 1, -3, -1, 0),
+        twoPhaseProperties.subDict
+        (
+            twoPhaseProperties.phase1Name()
+        ).lookup("k")
+    );
+
+    dimensionedScalar k2
+    (
+        "k",
+        dimensionSet(1, 1, -3, -1, 0),
+        twoPhaseProperties.subDict
+        (
+            twoPhaseProperties.phase2Name()
+        ).lookup("k")
+    );
+
+    dimensionedScalar Cv1
+    (
+        "Cv",
+        dimensionSet(0, 2, -2, -1, 0),
+        twoPhaseProperties.subDict
+        (
+            twoPhaseProperties.phase1Name()
+        ).lookup("Cv")
+    );
+
+    dimensionedScalar Cv2
+    (
+        "Cv",
+        dimensionSet(0, 2, -2, -1, 0),
+        twoPhaseProperties.subDict
+        (
+            twoPhaseProperties.phase2Name()
+        ).lookup("Cv")
+    );
+
+    autoPtr<phaseEquationOfState> eos1
+    (
+        phaseEquationOfState::New
+        (
+            twoPhaseProperties.subDict
+            (
+                twoPhaseProperties.phase1Name()
+            )
         )
     );
 
-    volScalarField rho1(rho10 + psi1*p);
-    volScalarField rho2(rho20 + psi2*p);
+    autoPtr<phaseEquationOfState> eos2
+    (
+        phaseEquationOfState::New
+        (
+            twoPhaseProperties.subDict
+            (
+                twoPhaseProperties.phase2Name()
+            )
+        )
+    );
+
+    volScalarField psi1
+    (
+        IOobject
+        (
+            "psi1",
+            runTime.timeName(),
+            mesh
+        ),
+        eos1->psi(p, T)
+    );
+    psi1.oldTime();
+
+    volScalarField psi2
+    (
+        IOobject
+        (
+            "psi2",
+            runTime.timeName(),
+            mesh
+        ),
+        eos2->psi(p, T)
+    );
+    psi2.oldTime();
+
+    dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    volScalarField rho1("rho1", eos1->rho(p, T));
+    volScalarField rho2("rho2", eos2->rho(p, T));
 
     volScalarField rho
     (

applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "wallHeatTransferFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    Tinf_(p.size(), 0.0),
+    alphaWall_(p.size(), 0.0)
+{
+    refValue() = 0.0;
+    refGrad() = 0.0;
+    valueFraction() = 0.0;
+}
+
+
+Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField
+(
+    const wallHeatTransferFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(ptf, p, iF, mapper),
+    Tinf_(ptf.Tinf_, mapper),
+    alphaWall_(ptf.alphaWall_, mapper)
+{}
+
+
+Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    Tinf_("Tinf", dict, p.size()),
+    alphaWall_("alphaWall", dict, p.size())
+{
+    refValue() = Tinf_;
+    refGrad() = 0.0;
+    valueFraction() = 0.0;
+
+    if (dict.found("value"))
+    {
+        fvPatchField<scalar>::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        evaluate();
+    }
+}
+
+
+Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField
+(
+    const wallHeatTransferFvPatchScalarField& tppsf
+)
+:
+    mixedFvPatchScalarField(tppsf),
+    Tinf_(tppsf.Tinf_),
+    alphaWall_(tppsf.alphaWall_)
+{}
+
+
+Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField
+(
+    const wallHeatTransferFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(tppsf, iF),
+    Tinf_(tppsf.Tinf_),
+    alphaWall_(tppsf.alphaWall_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::wallHeatTransferFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    scalarField::autoMap(m);
+    Tinf_.autoMap(m);
+    alphaWall_.autoMap(m);
+}
+
+
+void Foam::wallHeatTransferFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    mixedFvPatchScalarField::rmap(ptf, addr);
+
+    const wallHeatTransferFvPatchScalarField& tiptf =
+        refCast<const wallHeatTransferFvPatchScalarField>(ptf);
+
+    Tinf_.rmap(tiptf.Tinf_, addr);
+    alphaWall_.rmap(tiptf.alphaWall_, addr);
+}
+
+
+void Foam::wallHeatTransferFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const fvPatchScalarField& Cpw =
+        patch().lookupPatchField<volScalarField, scalar>("Cp");
+
+    const fvPatchScalarField& kByCpw =
+        patch().lookupPatchField<volScalarField, scalar>("kByCp");
+
+    valueFraction() =
+        1.0/
+        (
+            1.0
+          + Cpw*kByCpw*patch().deltaCoeffs()/alphaWall_
+        );
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::wallHeatTransferFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    Tinf_.writeEntry("Tinf", os);
+    alphaWall_.writeEntry("alphaWall", os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField(fvPatchScalarField, wallHeatTransferFvPatchScalarField);
+}
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.H

@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::wallHeatTransferFvPatchScalarField
+
+Description
+    Enthalpy boundary conditions for wall heat transfer
+
+SourceFiles
+    wallHeatTransferFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallHeatTransferFvPatchScalarField_H
+#define wallHeatTransferFvPatchScalarField_H
+
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class wallHeatTransferFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallHeatTransferFvPatchScalarField
+:
+    public mixedFvPatchScalarField
+{
+    // Private data
+
+        //- Tinf
+        scalarField Tinf_;
+
+        //- alphaWall
+        scalarField alphaWall_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("wallHeatTransfer");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        wallHeatTransferFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        wallHeatTransferFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given wallHeatTransferFvPatchScalarField
+        // onto a new patch
+        wallHeatTransferFvPatchScalarField
+        (
+            const wallHeatTransferFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        wallHeatTransferFvPatchScalarField
+        (
+            const wallHeatTransferFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new wallHeatTransferFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        wallHeatTransferFvPatchScalarField
+        (
+            const wallHeatTransferFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new wallHeatTransferFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return Tinf
+            const scalarField& Tinf() const
+            {
+                return Tinf_;
+            }
+
+            //- Return reference to Tinf to allow adjustment
+            scalarField& Tinf()
+            {
+                return Tinf_;
+            }
+
+            //- Return alphaWall
+            const scalarField& alphaWall() const
+            {
+                return alphaWall_;
+            }
+
+            //- Return reference to alphaWall to allow adjustment
+            scalarField& alphaWall()
+            {
+                return alphaWall_;
+            }
+
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -1,84 +1,97 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    rho1 = eos1->rho(p, T);
+    rho2 = eos2->rho(p, T);
 
-    tmp<fvScalarMatrix> p_rghEqnComp;
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    if (pimple.transonic())
+    tmp<fvScalarMatrix> p_rghEqnComp1;
+    tmp<fvScalarMatrix> p_rghEqnComp2;
+
+    //if (transonic)
+    //{
+    //}
+    //else
     {
-        p_rghEqnComp =
-        (
-            fvm::ddt(p_rgh)
-          + fvm::div(phi, p_rgh)
-          - fvm::Sp(fvc::div(phi), p_rgh)
-        );
-    }
-    else
-    {
-        p_rghEqnComp =
-        (
-            fvm::ddt(p_rgh)
-          + fvc::div(phi, p_rgh)
-          - fvc::Sp(fvc::div(phi), p_rgh)
-        );
+        surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
+        surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
+
+        p_rghEqnComp1 =
+            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+          + fvc::div(phid1, p_rgh)
+          - fvc::Sp(fvc::div(phid1), p_rgh);
+
+        p_rghEqnComp2 =
+            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+          + fvc::div(phid2, p_rgh)
+          - fvc::Sp(fvc::div(phid2), p_rgh);
     }
 
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    U = rAU*UEqn.H();
-
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    phi = phiU +
+    surfaceScalarField phig
+    (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    //thermo.rho() -= psi*p_rgh;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqnIncomp
         (
-            fvc::div(phi)
+            fvc::div(phiHbyA)
           - fvm::laplacian(rAUf, p_rgh)
         );
 
         solve
         (
             (
-                max(alpha1, scalar(0))*(psi1/rho1)
-              + max(alpha2, scalar(0))*(psi2/rho2)
+                (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
+              + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
             )
-           *p_rghEqnComp()
           + p_rghEqnIncomp,
             mesh.solver(p_rgh.select(pimple.finalInnerIter()))
         );
 
         if (pimple.finalNonOrthogonalIter())
         {
+            // Second part of thermodynamic density update
+            //thermo.rho() += psi*p_rgh;
+
             dgdt =
-                (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(p_rghEqnComp & p_rgh);
-            phi += p_rghEqnIncomp.flux();
-        }
-    }
-
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
-    p = max
             (
-        (p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
-       /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-        pMin
+                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
+              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
             );
 
-    rho1 = rho10 + psi1*p;
-    rho2 = rho20 + psi2*p;
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
+        }
+    }
+
+    p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+
+    rho1 = eos1->rho(p, T);
+    rho2 = eos2->rho(p, T);
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
     Info<< "min(p_rgh) " << min(p_rgh).value() << endl;

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/Make/files

@@ -0,0 +1,8 @@
+phaseEquationOfState/phaseEquationOfState.C
+phaseEquationOfState/newPhaseEquationOfState.C
+constant/constant.C
+linear/linear.C
+perfectFluid/perfectFluid.C
+adiabaticPerfectFluid/adiabaticPerfectFluid.C
+
+LIB = $(FOAM_LIBBIN)/libphaseEquationsOfState

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/Make/options

@@ -0,0 +1,6 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude
+
+LIB_LIBS = \
+    -lincompressibleTransportModels

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.C

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "adiabaticPerfectFluid.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace phaseEquationsOfState
+{
+    defineTypeNameAndDebug(adiabaticPerfectFluid, 0);
+
+    addToRunTimeSelectionTable
+    (
+        phaseEquationOfState,
+        adiabaticPerfectFluid,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseEquationsOfState::adiabaticPerfectFluid::adiabaticPerfectFluid
+(
+    const dictionary& dict
+)
+:
+    phaseEquationOfState(dict),
+    p0_("p0", dimPressure, dict.lookup("p0")),
+    rho0_("rho0", dimDensity, dict.lookup("rho0")),
+    gamma_("gamma", dimless, dict.lookup("gamma")),
+    B_("B", dimPressure, dict.lookup("B"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::phaseEquationsOfState::adiabaticPerfectFluid::~adiabaticPerfectFluid()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::phaseEquationsOfState::adiabaticPerfectFluid::rho
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rho",
+                p.time().timeName(),
+                p.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            rho0_*pow((p + B_)/(p0_ + B_), 1.0/gamma_)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::phaseEquationsOfState::adiabaticPerfectFluid::psi
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "psi",
+                p.time().timeName(),
+                p.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            (rho0_/(gamma_*(p0_ + B_)))
+           *pow((p + B_)/(p0_ + B_), 1.0/gamma_ - 1.0)
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phaseEquationsOfState::adiabaticPerfectFluid
+
+Description
+    AdiabaticPerfectFluid phase density model.
+
+SourceFiles
+    adiabaticPerfectFluid.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef adiabaticPerfectFluid_H
+#define adiabaticPerfectFluid_H
+
+#include "phaseEquationOfState.H"
+#include "dimensionedTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace phaseEquationsOfState
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class adiabaticPerfectFluid Declaration
+\*---------------------------------------------------------------------------*/
+
+class adiabaticPerfectFluid
+:
+    public phaseEquationOfState
+{
+    // Private data
+
+        //- Reference pressure
+        dimensionedScalar p0_;
+
+        //- Reference density
+        dimensionedScalar rho0_;
+
+        //- The isentropic exponent
+        dimensionedScalar gamma_;
+
+        //- Pressure offset for a stiffened gas
+        dimensionedScalar B_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("adiabaticPerfectFluid");
+
+
+    // Constructors
+
+        //- Construct from components
+        adiabaticPerfectFluid
+        (
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~adiabaticPerfectFluid();
+
+
+    // Member Functions
+
+        tmp<volScalarField> rho
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const;
+
+        tmp<volScalarField> psi
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace phaseEquationsOfState
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/constant/constant.C

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "constant.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace phaseEquationsOfState
+{
+    defineTypeNameAndDebug(constant, 0);
+
+    addToRunTimeSelectionTable
+    (
+        phaseEquationOfState,
+        constant,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseEquationsOfState::constant::constant
+(
+    const dictionary& dict
+)
+:
+    phaseEquationOfState(dict),
+    rho_("rho", dimDensity, dict.lookup("rho"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::phaseEquationsOfState::constant::~constant()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::phaseEquationsOfState::constant::rho
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rho",
+                p.time().timeName(),
+                p.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            p.mesh(),
+            rho_
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::phaseEquationsOfState::constant::psi
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "psi",
+                p.time().timeName(),
+                p.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            p.mesh(),
+            dimensionedScalar("psi", dimDensity/dimPressure, 0)
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/constant/constant.H

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phaseEquationsOfState::constant
+
+Description
+    Constant phase density model.
+
+SourceFiles
+    constant.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constant_H
+#define constant_H
+
+#include "phaseEquationOfState.H"
+#include "dimensionedTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace phaseEquationsOfState
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class constant Declaration
+\*---------------------------------------------------------------------------*/
+
+class constant
+:
+    public phaseEquationOfState
+{
+    // Private data
+
+        //- The constant density of the phase
+        dimensionedScalar rho_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constant");
+
+
+    // Constructors
+
+        //- Construct from components
+        constant
+        (
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~constant();
+
+
+    // Member Functions
+
+        tmp<volScalarField> rho
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const;
+
+        tmp<volScalarField> psi
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace phaseEquationsOfState
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/linear/linear.C

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "linear.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace phaseEquationsOfState
+{
+    defineTypeNameAndDebug(linear, 0);
+
+    addToRunTimeSelectionTable
+    (
+        phaseEquationOfState,
+        linear,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseEquationsOfState::linear::linear
+(
+    const dictionary& dict
+)
+:
+    phaseEquationOfState(dict),
+    rho0_("rho0", dimDensity, dict.lookup("rho0")),
+    psi_("psi", dimDensity/dimPressure, dict.lookup("psi"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::phaseEquationsOfState::linear::~linear()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::phaseEquationsOfState::linear::rho
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rho",
+                p.time().timeName(),
+                p.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            rho0_ + psi_*p
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::phaseEquationsOfState::linear::psi
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "psi",
+                p.time().timeName(),
+                p.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            p.mesh(),
+            psi_
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/linear/linear.H

@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phaseEquationsOfState::linear
+
+Description
+    Linear phase density model.
+
+SourceFiles
+    linear.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linear_H
+#define linear_H
+
+#include "phaseEquationOfState.H"
+#include "dimensionedTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace phaseEquationsOfState
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class linear Declaration
+\*---------------------------------------------------------------------------*/
+
+class linear
+:
+    public phaseEquationOfState
+{
+    // Private data
+
+        //- The reference density of the phase
+        dimensionedScalar rho0_;
+
+        //- The constant compressibility of the phase
+        dimensionedScalar psi_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("linear");
+
+
+    // Constructors
+
+        //- Construct from components
+        linear
+        (
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~linear();
+
+
+    // Member Functions
+
+        tmp<volScalarField> rho
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const;
+
+        tmp<volScalarField> psi
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace phaseEquationsOfState
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/perfectFluid/perfectFluid.C

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "perfectFluid.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace phaseEquationsOfState
+{
+    defineTypeNameAndDebug(perfectFluid, 0);
+
+    addToRunTimeSelectionTable
+    (
+        phaseEquationOfState,
+        perfectFluid,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseEquationsOfState::perfectFluid::perfectFluid
+(
+    const dictionary& dict
+)
+:
+    phaseEquationOfState(dict),
+    rho0_("rho0", dimDensity, dict.lookup("rho0")),
+    R_("R", dimensionSet(0, 2, -2, -1, 0), dict.lookup("R"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::phaseEquationsOfState::perfectFluid::~perfectFluid()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::phaseEquationsOfState::perfectFluid::rho
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rho",
+                p.time().timeName(),
+                p.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            rho0_ + psi(p, T)*p
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::phaseEquationsOfState::perfectFluid::psi
+(
+    const volScalarField& p,
+    const volScalarField& T
+) const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "psi",
+                p.time().timeName(),
+                p.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            1.0/(R_*T)
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/perfectFluid/perfectFluid.H

@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phaseEquationsOfState::perfectFluid
+
+Description
+    PerfectFluid phase density model.
+
+SourceFiles
+    perfectFluid.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef perfectFluid_H
+#define perfectFluid_H
+
+#include "phaseEquationOfState.H"
+#include "dimensionedTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace phaseEquationsOfState
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class perfectFluid Declaration
+\*---------------------------------------------------------------------------*/
+
+class perfectFluid
+:
+    public phaseEquationOfState
+{
+    // Private data
+
+        //- The reference density of the phase
+        dimensionedScalar rho0_;
+
+        //- The fluid constant of the phase
+        dimensionedScalar R_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("perfectFluid");
+
+
+    // Constructors
+
+        //- Construct from components
+        perfectFluid
+        (
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~perfectFluid();
+
+
+    // Member Functions
+
+        tmp<volScalarField> rho
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const;
+
+        tmp<volScalarField> psi
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace phaseEquationsOfState
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/newPhaseEquationOfState.C

@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phaseEquationOfState.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::phaseEquationOfState> Foam::phaseEquationOfState::New
+(
+    const dictionary& dict
+)
+{
+    word phaseEquationOfStateType
+    (
+        dict.subDict("equationOfState").lookup("type")
+    );
+
+    Info<< "Selecting phaseEquationOfState "
+        << phaseEquationOfStateType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(phaseEquationOfStateType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn("phaseEquationOfState::New")
+           << "Unknown phaseEquationOfStateType type "
+           << phaseEquationOfStateType << endl << endl
+           << "Valid phaseEquationOfState types are : " << endl
+           << dictionaryConstructorTablePtr_->sortedToc()
+           << exit(FatalError);
+    }
+
+    return cstrIter()(dict.subDict("equationOfState"));
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/phaseEquationOfState.C

@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phaseEquationOfState.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(phaseEquationOfState, 0);
+    defineRunTimeSelectionTable(phaseEquationOfState, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseEquationOfState::phaseEquationOfState
+(
+    const dictionary& dict
+)
+:
+    dict_(dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::phaseEquationOfState::~phaseEquationOfState()
+{}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/phaseEquationOfState.H

@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phaseEquationOfState
+
+Description
+    A2stract base-class for dispersed-phase particle diameter models.
+
+SourceFiles
+    phaseEquationOfState.C
+    newDiameterModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseEquationOfState_H
+#define phaseEquationOfState_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "dictionary.H"
+#include "volFieldsFwd.H"
+#include "runTimeSelectionTables.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class phaseEquationOfState Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseEquationOfState
+{
+protected:
+
+    // Protected data
+
+        const dictionary& dict_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("phaseEquationOfState");
+
+
+    // Declare runtime construction
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            phaseEquationOfState,
+            dictionary,
+            (
+                const dictionary& dict
+            ),
+            (dict)
+        );
+
+
+    // Constructors
+
+        phaseEquationOfState
+        (
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~phaseEquationOfState();
+
+
+    // Selectors
+
+        static autoPtr<phaseEquationOfState> New
+        (
+            const dictionary& dict
+        );
+
+
+    // Member Functions
+
+        //- Return the phase density
+        virtual tmp<volScalarField> rho
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const = 0;
+
+        //- Return the phase compressibility
+        virtual tmp<volScalarField> psi
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H

@@ -1,29 +1,35 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
-    surfaceScalarField phiU
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiU);
+    mrfZones.relativeFlux(phiHbyA);
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi = phiU +
+    surfaceScalarField phig
+    (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -32,12 +38,12 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
-        }
-    }
+            phi = phiHbyA - p_rghEqn.flux();
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
+        }
+    }
 
     #include "continuityErrs.H"
 

applications/solvers/multiphase/interFoam/createFields.H

@@ -12,20 +12,6 @@
         mesh
     );
 
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -46,6 +32,8 @@
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi);
 
+    volScalarField& alpha1(twoPhaseProperties.alpha1());
+
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
 

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -1,31 +1,39 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    phiAbs =
-        (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phiAbs);
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+    );
 
     if (p_rgh.needReference())
     {
-        fvc::makeRelative(phiAbs, U);
-        adjustPhi(phiAbs, U, p_rgh);
-        fvc::makeAbsolute(phiAbs, U);
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p_rgh);
+        fvc::makeAbsolute(phiHbyA, U);
     }
 
-    phi = phiAbs +
+    surfaceScalarField phig
+    (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -34,12 +42,12 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
-        }
-    }
+            phi = phiHbyA - p_rghEqn.flux();
 
-    U += rAU*fvc::reconstruct((phi - phiAbs)/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
+        }
+    }
 
     #include "continuityErrs.H"
 

applications/solvers/multiphase/interFoam/pEqn.H

@@ -1,28 +1,34 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
-    surfaceScalarField phiU
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi = phiU +
+    surfaceScalarField phig
+    (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -31,12 +37,12 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
-        }
-    }
+            phi = phiHbyA - p_rghEqn.flux();
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
+        }
+    }
 
     #include "continuityErrs.H"
 

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -1,24 +1,28 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi =
-        phiU
-      + (
+    surfaceScalarField phig
+    (
+        (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
@@ -28,7 +32,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvc::div(phi) - fvm::laplacian(rAUf, p_rgh)
+            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
           - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
         );
 
@@ -38,12 +42,12 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += p_rghEqn.flux();
-        }
-    }
+            phi = phiHbyA + p_rghEqn.flux();
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
+        }
+    }
 
     #include "continuityErrs.H"
 

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H

@@ -1,31 +1,32 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
-      //+ fvc::ddtPhiCorr(rAU, rho, U, phi)
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiU);
+    mrfZones.relativeFlux(phiHbyA);
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi =
-        phiU
-      + (
-            mixture.surfaceTensionForce()
-          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -34,12 +35,12 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
-        }
-    }
+            phi = phiHbyA - p_rghEqn.flux();
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
+        }
+    }
 
     #include "continuityErrs.H"
 

applications/solvers/multiphase/multiphaseInterFoam/pEqn.H

@@ -1,29 +1,34 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi = phiU +
+    surfaceScalarField phig
+    (
         (
             mixture.surfaceTensionForce()
           - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -32,12 +37,12 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
-        }
-    }
+            phi = phiHbyA - p_rghEqn.flux();
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
+        }
+    }
 
     #include "continuityErrs.H"
 

applications/solvers/multiphase/settlingFoam/pEqn.H

@@ -1,27 +1,33 @@
-volScalarField rAU(1.0/UEqn.A());
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
 
-surfaceScalarField rAUf
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rho*rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-    "(rho*(1|A(U)))",
-    fvc::interpolate(rho)*fvc::interpolate(rAU)
-);
-
-U = rAU*UEqn.H();
-phi =
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-       (fvc::interpolate(U) & mesh.Sf())
+           (fvc::interpolate(HbyA) & mesh.Sf())
          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        )
     );
 
-surfaceScalarField phiU("phiU", phi);
-phi -= ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-        fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -30,7 +36,10 @@ while (pimple.correctNonOrthogonal())
 
         if (pimple.finalNonOrthogonalIter())
         {
-        phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
@@ -49,6 +58,4 @@ if (p_rgh.needReference())
 
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
-
-U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-U.correctBoundaryConditions();
+}

applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H

@@ -12,20 +12,6 @@
         mesh
     );
 
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -45,6 +31,8 @@
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi);
 
+    volScalarField& alpha1(twoPhaseProperties.alpha1());
+
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
 

applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H

@@ -1,24 +1,31 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
-    surfaceScalarField phiU
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi = phiU - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -27,12 +34,12 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
-        }
-    }
+            phi = phiHbyA - p_rghEqn.flux();
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
+        }
+    }
 
     #include "continuityErrs.H"
 

applications/test/LduMatrix/LduMatrixTest.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "LduMatrix.H"
+#include "diagTensorField.H"
+#include "TPCG.H"
+#include "TPBiCG.H"
+#include "NoPreconditioner.H"
+
+using namespace Foam;
+
+typedef Foam::LduMatrix<vector, diagTensor, scalar>
+    lduVectorMatrix;
+defineNamedTemplateTypeNameAndDebug(lduVectorMatrix, 0);
+
+typedef Foam::DiagonalSolver<vector, diagTensor, scalar>
+    lduVectorDiagonalSolver;
+defineNamedTemplateTypeNameAndDebug(lduVectorDiagonalSolver, 0);
+
+template<>
+const vector lduVectorMatrix::great_(1e15, 1e15, 1e15);
+
+template<>
+const vector lduVectorMatrix::small_(1e-15, 1e-15, 1e-15);
+
+namespace Foam
+{
+    typedef LduMatrix<vector, diagTensor, scalar>::preconditioner
+        lduVectorPreconditioner;
+    defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, asymMatrix);
+
+    typedef LduMatrix<vector, diagTensor, scalar>::smoother
+        lduVectorSmoother;
+    defineTemplateRunTimeSelectionTable(lduVectorSmoother, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorSmoother, asymMatrix);
+
+    typedef LduMatrix<vector, diagTensor, scalar>::solver
+        lduVectorSolver;
+    defineTemplateRunTimeSelectionTable(lduVectorSolver, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorSolver, asymMatrix);
+
+    typedef TPCG<vector, diagTensor, scalar> TPCGVector;
+    defineNamedTemplateTypeNameAndDebug(TPCGVector, 0);
+
+    LduMatrix<vector, diagTensor, scalar>::solver::
+        addsymMatrixConstructorToTable<TPCGVector>
+        addTPCGSymMatrixConstructorToTable_;
+
+    typedef TPBiCG<vector, diagTensor, scalar> TPBiCGVector;
+    defineNamedTemplateTypeNameAndDebug(TPBiCGVector, 0);
+
+    LduMatrix<vector, diagTensor, scalar>::solver::
+        addasymMatrixConstructorToTable<TPBiCGVector>
+        addTPBiCGSymMatrixConstructorToTable_;
+
+    typedef NoPreconditioner<vector, diagTensor, scalar> NoPreconditionerVector;
+    defineNamedTemplateTypeNameAndDebug(NoPreconditionerVector, 0);
+
+    LduMatrix<vector, diagTensor, scalar>::preconditioner::
+        addsymMatrixConstructorToTable<NoPreconditionerVector>
+        addNoPreconditionerSymMatrixConstructorToTable_;
+
+    LduMatrix<vector, diagTensor, scalar>::preconditioner::
+        addasymMatrixConstructorToTable<NoPreconditionerVector>
+        addNoPreconditionerAsymMatrixConstructorToTable_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    volVectorField psi
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    lduVectorMatrix testMatrix(mesh);
+    testMatrix.diag() = 2*pTraits<diagTensor>::one;
+    testMatrix.source() = pTraits<vector>::one;
+    testMatrix.upper() = 0.1;
+    testMatrix.lower() = -0.1;
+
+    Info<< testMatrix << endl;
+
+    FieldField<Field, scalar> boundaryCoeffs(0);
+    FieldField<Field, scalar> internalCoeffs(0);
+
+    autoPtr<lduVectorMatrix::solver> testMatrixSolver =
+    lduVectorMatrix::solver::New
+    (
+        psi.name(),
+        testMatrix,
+        boundaryCoeffs,
+        internalCoeffs,
+        psi.boundaryField().interfaces(),
+        IStringStream
+        (
+            "PBiCG"
+            "{"
+            "    preconditioner   none;"
+            "    tolerance        (1e-05 1e-05 1e-05);"
+            "    relTol           (0 0 0);"
+            "}"
+        )()
+    );
+
+    lduVectorMatrix::solverPerformance solverPerf =
+        testMatrixSolver->solve(psi);
+
+    solverPerf.print();
+
+    Info<< psi << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/LduMatrix/LduMatrixTest2.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "LduMatrix.H"
+#include "tensorField.H"
+#include "TPCG.H"
+#include "TPBiCG.H"
+#include "NoPreconditioner.H"
+
+using namespace Foam;
+
+typedef Foam::LduMatrix<vector, tensor, scalar>
+    lduVectorMatrix;
+defineNamedTemplateTypeNameAndDebug(lduVectorMatrix, 0);
+
+typedef Foam::DiagonalSolver<vector, tensor, scalar>
+    lduVectorDiagonalSolver;
+defineNamedTemplateTypeNameAndDebug(lduVectorDiagonalSolver, 0);
+
+template<>
+const vector lduVectorMatrix::great_(1e15, 1e15, 1e15);
+
+template<>
+const vector lduVectorMatrix::small_(1e-15, 1e-15, 1e-15);
+
+namespace Foam
+{
+    typedef LduMatrix<vector, tensor, scalar>::preconditioner
+        lduVectorPreconditioner;
+    defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, asymMatrix);
+
+    typedef LduMatrix<vector, tensor, scalar>::smoother
+        lduVectorSmoother;
+    defineTemplateRunTimeSelectionTable(lduVectorSmoother, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorSmoother, asymMatrix);
+
+    typedef LduMatrix<vector, tensor, scalar>::solver
+        lduVectorSolver;
+    defineTemplateRunTimeSelectionTable(lduVectorSolver, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorSolver, asymMatrix);
+
+    typedef TPCG<vector, tensor, scalar> TPCGVector;
+    defineNamedTemplateTypeNameAndDebug(TPCGVector, 0);
+
+    LduMatrix<vector, tensor, scalar>::solver::
+        addsymMatrixConstructorToTable<TPCGVector>
+        addTPCGSymMatrixConstructorToTable_;
+
+    typedef TPBiCG<vector, tensor, scalar> TPBiCGVector;
+    defineNamedTemplateTypeNameAndDebug(TPBiCGVector, 0);
+
+    LduMatrix<vector, tensor, scalar>::solver::
+        addasymMatrixConstructorToTable<TPBiCGVector>
+        addTPBiCGSymMatrixConstructorToTable_;
+
+    typedef NoPreconditioner<vector, tensor, scalar> NoPreconditionerVector;
+    defineNamedTemplateTypeNameAndDebug(NoPreconditionerVector, 0);
+
+    LduMatrix<vector, tensor, scalar>::preconditioner::
+        addsymMatrixConstructorToTable<NoPreconditionerVector>
+        addNoPreconditionerSymMatrixConstructorToTable_;
+
+    LduMatrix<vector, tensor, scalar>::preconditioner::
+        addasymMatrixConstructorToTable<NoPreconditionerVector>
+        addNoPreconditionerAsymMatrixConstructorToTable_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    volVectorField psi
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    lduVectorMatrix testMatrix(mesh);
+    testMatrix.diag() = 2*I;
+    testMatrix.source() = pTraits<vector>::one;
+    testMatrix.upper() = 0.1;
+    testMatrix.lower() = -0.1;
+
+    Info<< testMatrix << endl;
+
+    FieldField<Field, scalar> boundaryCoeffs(0);
+    FieldField<Field, scalar> internalCoeffs(0);
+
+    autoPtr<lduVectorMatrix::solver> testMatrixSolver =
+    lduVectorMatrix::solver::New
+    (
+        psi.name(),
+        testMatrix,
+        //boundaryCoeffs,
+        //internalCoeffs,
+        //psi.boundaryField().interfaces(),
+        IStringStream
+        (
+            "PBiCG"
+            "{"
+            "    preconditioner   none;"
+            "    tolerance        (1e-05 1e-05 1e-05);"
+            "    relTol           (0 0 0);"
+            "}"
+        )()
+    );
+
+    lduVectorMatrix::solverPerformance solverPerf =
+        testMatrixSolver->solve(psi);
+
+    solverPerf.print();
+
+    Info<< psi << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/LduMatrix/LduMatrixTest3.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+Application
+    icoFoam
+
+Description
+    Transient solver for incompressible, laminar flow of Newtonian fluids.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "LduMatrix.H"
+#include "diagTensorField.H"
+
+typedef LduMatrix<vector, scalar, scalar> lduVectorMatrix;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readPISOControls.H"
+        #include "CourantNo.H"
+
+        fvVectorMatrix UEqn
+        (
+            fvm::ddt(U)
+          + fvm::div(phi, U)
+          - fvm::laplacian(nu, U)
+        );
+
+        fvVectorMatrix UEqnp(UEqn == -fvc::grad(p));
+
+        lduVectorMatrix U3Eqnp(mesh);
+        U3Eqnp.diag() = UEqnp.diag();
+        U3Eqnp.upper() = UEqnp.upper();
+        U3Eqnp.lower() = UEqnp.lower();
+        U3Eqnp.source() = UEqnp.source();
+
+        UEqnp.addBoundaryDiag(U3Eqnp.diag(), 0);
+        UEqnp.addBoundarySource(U3Eqnp.source(), false);
+
+        autoPtr<lduVectorMatrix::solver> U3EqnpSolver =
+            lduVectorMatrix::solver::New
+            (
+                U.name(),
+                U3Eqnp,
+                dictionary
+                (
+                    IStringStream
+                    (
+                        "{"
+                        "    solver           PBiCG;"
+                        "    preconditioner   DILU;"
+                        "    tolerance        (1e-5 1e-5 1);"
+                        "    relTol           (0 0 0);"
+                        "}"
+                    )()
+                )
+            );
+
+        U3EqnpSolver->solve(U).print(Info);
+
+        // --- PISO loop
+
+        for (int corr=0; corr<nCorr; corr++)
+        {
+            volScalarField rAU = 1.0/UEqn.A();
+
+            U = rAU*UEqn.H();
+            phi = (fvc::interpolate(U) & mesh.Sf())
+                + fvc::ddtPhiCorr(rAU, U, phi);
+
+            adjustPhi(phi, U, p);
+
+            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            {
+                fvScalarMatrix pEqn
+                (
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
+                );
+
+                pEqn.setReference(pRefCell, pRefValue);
+                pEqn.solve();
+
+                if (nonOrth == nNonOrthCorr)
+                {
+                    phi -= pEqn.flux();
+                }
+            }
+
+            #include "continuityErrs.H"
+
+            U -= rAU*fvc::grad(p);
+            U.correctBoundaryConditions();
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/LduMatrix/Make/files

@@ -0,0 +1,3 @@
+LduMatrixTest3.C
+
+EXE = $(FOAM_USER_APPBIN)/LduMatrixTest

applications/test/LduMatrix/Make/options

@@ -0,0 +1,4 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = -lfiniteVolume

applications/test/LduMatrix/createFields.H

@@ -0,0 +1,55 @@
+    Info<< "Reading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    dimensionedScalar nu
+    (
+        transportProperties.lookup("nu")
+    );
+
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+#   include "createPhi.H"
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

applications/test/PisoFoam/Make/files

@@ -0,0 +1,3 @@
+PisoFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/PisoFoam

applications/test/PisoFoam/Make/options

@@ -0,0 +1,13 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools

applications/test/PisoFoam/PisoFoam.C

@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+Application
+    pisoFoam
+
+Description
+    Transient solver for incompressible flow.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
+
+#include "LduMatrix.H"
+#include "diagTensorField.H"
+
+typedef LduMatrix<vector, scalar, scalar> lduVectorMatrix;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readPISOControls.H"
+        #include "CourantNo.H"
+
+        // Pressure-velocity PISO corrector
+        {
+            // Momentum predictor
+
+            fvVectorMatrix UEqn
+            (
+                fvm::ddt(U)
+              + fvm::div(phi, U)
+              + turbulence->divDevReff(U)
+            );
+
+            //UEqn.relax();
+
+            fvVectorMatrix UEqnp(UEqn == -fvc::grad(p));
+
+            lduVectorMatrix U3Eqnp(mesh);
+            U3Eqnp.diag() = UEqnp.diag();
+            U3Eqnp.upper() = UEqnp.upper();
+            U3Eqnp.lower() = UEqnp.lower();
+            U3Eqnp.source() = UEqnp.source();
+
+            UEqnp.addBoundaryDiag(U3Eqnp.diag(), 0);
+            UEqnp.addBoundarySource(U3Eqnp.source(), false);
+
+            U3Eqnp.interfaces() = U.boundaryField().interfaces();
+            U3Eqnp.interfacesUpper() = UEqnp.boundaryCoeffs().component(0);
+            U3Eqnp.interfacesLower() = UEqnp.internalCoeffs().component(0);
+
+            autoPtr<lduVectorMatrix::solver> U3EqnpSolver =
+            lduVectorMatrix::solver::New
+            (
+                U.name(),
+                U3Eqnp,
+                dictionary
+                (
+                    IStringStream
+                    (
+                        "{"
+                        "    /*solver          SmoothSolver;*/"
+                        "    smoother        GaussSeidel;"
+                        "    solver           PBiCCCG;"
+                        "    preconditioner   none;"
+                        "    tolerance        (1e-7 1e-7 1);"
+                        "    relTol           (0 0 0);"
+                        "}"
+                    )()
+                )
+            );
+
+            //for (int i=0; i<3; i++)
+            {
+                U3EqnpSolver->solve(U).print(Info);
+                U.correctBoundaryConditions();
+            }
+            //solve(UEqnp);
+
+            // --- PISO loop
+
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                volScalarField rAU(1.0/UEqn.A());
+
+                volVectorField HbyA("HbyA", U);
+                HbyA = rAU*UEqn.H();
+                surfaceScalarField phiHbyA
+                (
+                    "phiHbyA",
+                    (fvc::interpolate(HbyA) & mesh.Sf())
+                  + fvc::ddtPhiCorr(rAU, U, phi)
+                );
+
+                adjustPhi(phiHbyA, U, p);
+
+                // Non-orthogonal pressure corrector loop
+                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                {
+                    // Pressure corrector
+
+                    fvScalarMatrix pEqn
+                    (
+                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+                    );
+
+                    pEqn.setReference(pRefCell, pRefValue);
+
+                    if
+                    (
+                        corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
+                    )
+                    {
+                        pEqn.solve(mesh.solver("pFinal"));
+                    }
+                    else
+                    {
+                        pEqn.solve();
+                    }
+
+                    if (nonOrth == nNonOrthCorr)
+                    {
+                        phi = phiHbyA - pEqn.flux();
+                    }
+                }
+
+                #include "continuityErrs.H"
+
+                U = HbyA - rAU*fvc::grad(p);
+                U.correctBoundaryConditions();
+            }
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/PisoFoam/createFields.H

@@ -0,0 +1,42 @@
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );

applications/test/extendedStencil/Test-ExtendedStencil2.C

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 <http://www.gnu.org/licenses/>.
+
+Application
+    testExtendedStencil2
+
+Description
+    Test app for determining extended cell-to-cell stencil.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "Time.H"
+#include "OFstream.H"
+#include "meshTools.H"
+
+#include "CFCCellToCellStencil.H"
+
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void writeStencilOBJ
+(
+    const fileName& fName,
+    const point& fc,
+    const List<point>& stencilCc
+)
+{
+    OFstream str(fName);
+    label vertI = 0;
+
+    meshTools::writeOBJ(str, fc);
+    vertI++;
+
+    forAll(stencilCc, i)
+    {
+        meshTools::writeOBJ(str, stencilCc[i]);
+        vertI++;
+        str << "l 1 " << vertI << nl;
+    }
+}
+
+
+// Stats
+void writeStencilStats(const labelListList& stencil)
+{
+    label sumSize = 0;
+    label nSum = 0;
+    label minSize = labelMax;
+    label maxSize = labelMin;
+
+    forAll(stencil, i)
+    {
+        const labelList& sCells = stencil[i];
+
+        if (sCells.size() > 0)
+        {
+            sumSize += sCells.size();
+            nSum++;
+            minSize = min(minSize, sCells.size());
+            maxSize = max(maxSize, sCells.size());
+        }
+    }
+    reduce(sumSize, sumOp<label>());
+    reduce(nSum, sumOp<label>());
+    sumSize /= nSum;
+
+    reduce(minSize, minOp<label>());
+    reduce(maxSize, maxOp<label>());
+
+    Info<< "Stencil size :" << nl
+        << "    average : " << sumSize << nl
+        << "    min     : " << minSize << nl
+        << "    max     : " << maxSize << nl
+        << endl;
+}
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "addTimeOptions.H"
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    // Get times list
+    instantList Times = runTime.times();
+#   include "checkTimeOptions.H"
+    runTime.setTime(Times[startTime], startTime);
+#   include "createMesh.H"
+
+
+    //---- CELL CENTRED STENCIL -----
+
+    // Centred, cell stencil
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        // Full stencil. This is per local cell a set of global indices, either
+        // into cells or also boundary faces.
+        CFCCellToCellStencil stencil(mesh);
+
+        // Construct exchange map and renumber stencil
+        List<Map<label> > compactMap(Pstream::nProcs());
+        mapDistribute map
+        (
+            stencil.globalNumbering(),
+            stencil,
+            compactMap
+        );
+
+
+        // Print some stats
+        Info<< "cellFaceCell:" << endl;
+        writeStencilStats(stencil);
+
+
+        // Collect the data in stencil form
+        List<List<vector> > stencilPoints;
+        {
+            const volVectorField& fld = mesh.C();
+
+            // 1. Construct cell data in compact addressing
+            List<point> compactFld(map.constructSize(), pTraits<point>::zero);
+
+            // Insert my internal values
+            forAll(fld, cellI)
+            {
+                compactFld[cellI] = fld[cellI];
+            }
+            // Insert my boundary values
+            label nCompact = fld.size();
+            forAll(fld.boundaryField(), patchI)
+            {
+                const fvPatchField<vector>& pfld = fld.boundaryField()[patchI];
+
+                forAll(pfld, i)
+                {
+                    compactFld[nCompact++] = pfld[i];
+                }
+            }
+
+            // Do all swapping
+            map.distribute(compactFld);
+
+            // 2. Pull to stencil
+            stencilPoints.setSize(stencil.size());
+
+            forAll(stencil, cellI)
+            {
+                const labelList& compactCells = stencil[cellI];
+
+                stencilPoints[cellI].setSize(compactCells.size());
+
+                forAll(compactCells, i)
+                {
+                    stencilPoints[cellI][i] = compactFld[compactCells[i]];
+                }
+            }
+        }
+
+
+        forAll(stencilPoints, cellI)
+        {
+            writeStencilOBJ
+            (
+                runTime.path()/"centredCell" + Foam::name(cellI) + ".obj",
+                mesh.cellCentres()[cellI],
+                stencilPoints[cellI]
+            );
+        }
+    }
+
+    Pout<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/syncTools/Test-syncTools.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
@@ -214,7 +214,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
 
         forAll(localPoints, i)
         {
-            const point pt = localPoints[i] + 1E-4*rndGen.vector01();
+            const point pt = localPoints[i] + 1e-4*rndGen.vector01();
 
             label meshPointI = allBoundary.meshPoints()[i];
 
@@ -299,7 +299,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         {
             const edge& e = edges[i];
 
-            const point pt = e.centre(localPoints) + 1E-4*rndGen.vector01();
+            const point pt = e.centre(localPoints) + 1e-4*rndGen.vector01();
 
             label meshEdgeI = meshEdges[i];
 

applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.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
@@ -58,7 +58,7 @@ using namespace Foam;
 
 
 // Max cos angle for edges to be considered aligned with axis.
-static const scalar edgeTol = 1E-3;
+static const scalar edgeTol = 1e-3;
 
 
 void writeSet(const cellSet& cells, const string& msg)

applications/utilities/mesh/advanced/collapseEdges/collapseEdges.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
@@ -528,7 +528,7 @@ int main(int argc, char *argv[])
                 (
                     mesh,
                     boundaryPoint,
-                    1E-9,       // factor of largest face area
+                    1e-9,       // factor of largest face area
                     5,          // factor between smallest and largest edge on
                                 // face
                     collapser

applications/utilities/mesh/advanced/selectCells/edgeStats.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
@@ -31,7 +31,7 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::scalar Foam::edgeStats::edgeTol_ = 1E-3;
+const Foam::scalar Foam::edgeStats::edgeTol_ = 1e-3;
 
 
 

applications/utilities/mesh/generation/cv2DMesh/Make/options

@@ -10,7 +10,6 @@ EXE_INC = \
     ${EXE_NDEBUG} \
     ${CGAL_INC} \
     -I$(FOAM_APP)/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/lnInclude \
-    -I$(FOAM_APP)/utilities/mesh/generation/extrude/extrudeModel/lnInclude \
     -IconformalVoronoi2DMesh/lnInclude \
     -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -18,10 +17,13 @@ EXE_INC = \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/mesh/extrudeModel/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude
 
 EXE_LIBS = \
     $(CGAL_LIBS) \
+    -lboost_thread \
+    -lmpfr \
     -lextrude2DMesh \
     -lextrudeModel \
     -lcv2DMesh \

applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.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
@@ -209,7 +209,7 @@ int main(int argc, char *argv[])
             const pointField& points = pMesh.points();
 
             const boundBox& bb = pMesh.bounds();
-            const scalar mergeDim = 1E-4 * bb.minDim();
+            const scalar mergeDim = 1e-4 * bb.minDim();
 
             forAll(edges, edgeI)
             {

applications/utilities/mesh/generation/cvMesh/Make/options

@@ -22,6 +22,8 @@ EXE_INC = \
 
 EXE_LIBS = \
     $(CGAL_LIBS) \
+    -lboost_thread \
+    -lmpfr \
     -lconformalVoronoiMesh \
     -lmeshTools \
     -ldecompositionMethods \

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -45,8 +45,6 @@ SourceFiles
 #ifndef conformalVoronoiMesh_H
 #define conformalVoronoiMesh_H
 
-#define CGAL_INEXACT
-
 #include "CGALTriangulation3Ddefs.H"
 #include "uint.H"
 #include "ulong.H"

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -561,7 +561,7 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
     scalar closenessTolerance = cvMeshControls().mergeClosenessCoeff();
 
     // Absolute distance for points to be considered coincident. Bit adhoc
-    // but points were seen with distSqr ~ 1E-30 which is SMALL^2. Add a few
+    // but points were seen with distSqr ~ 1e-30 which is SMALL^2. Add a few
     // digits to account for truncation errors.
     scalar coincidentDistanceSqr = sqr
     (

applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/Make/options

@@ -19,6 +19,8 @@ EXE_INC = \
 
 EXE_LIBS = \
     $(CGAL_LIBS) \
+    -lboost_thread \
+    -lmpfr \
     -lconformalVoronoiMesh \
     -ldecompositionMethods /* -L$(FOAM_LIBBIN)/dummy -lscotchDecomp */ \
     -ledgeMesh \

applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/cvMeshBackgroundMesh.C

@@ -51,7 +51,7 @@ using namespace Foam;
 
 // Tolerance (as fraction of the bounding box). Needs to be fairly lax since
 // usually meshes get written with limited precision (6 digits)
-static const scalar defaultMergeTol = 1E-6;
+static const scalar defaultMergeTol = 1e-6;
 
 // Get merging distance when matching face centres
 scalar getMergeDistance
@@ -397,7 +397,7 @@ int main(int argc, char *argv[])
         "mergeTol",
         "scalar",
         "specify the merge distance relative to the bounding box size "
-        "(default 1E-6)"
+        "(default 1e-6)"
     );
 
     #include "setRootCase.H"

applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Make/options

@@ -14,6 +14,8 @@ EXE_INC = \
 
 EXE_LIBS = \
     $(CGAL_LIBS) \
+    -lboost_thread \
+    -lmpfr \
     -L$(FASTDUALOCTREE_SRC_PATH) -lperf_main \
     -lGL \
     -lconformalVoronoiMesh \

applications/utilities/mesh/generation/extrude/Allwclean

@@ -2,7 +2,6 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-wclean libso extrudeModel
 wclean extrudeMesh
 wclean extrudeToRegionMesh
 

applications/utilities/mesh/generation/extrude/Allwmake

@@ -2,10 +2,7 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-wmake libso extrudeModel
 wmake extrudeMesh
-
-wmake libso extrudeToRegionMesh/createShellMesh
 wmake extrudeToRegionMesh
 
 

applications/utilities/mesh/generation/extrude/extrudeMesh/Make/options

@@ -1,10 +1,10 @@
 EXE_INC = \
     -IextrudedMesh \
-    -I../extrudeModel/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/mesh/extrudeModel/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.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
@@ -666,7 +666,7 @@ int main(int argc, char *argv[])
 
     const boundBox& bb = mesh.bounds();
     const vector span = bb.span();
-    const scalar mergeDim = 1E-4 * bb.minDim();
+    const scalar mergeDim = 1e-4 * bb.minDim();
 
     Info<< "Mesh bounding box : " << bb << nl
         << "        with span : " << span << nl

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/options

@@ -1,13 +1,11 @@
 EXE_INC = \
-    -I../extrudeModel/lnInclude \
-    -IcreateShellMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/mesh/extrudeModel/lnInclude
 
 EXE_LIBS = \
-    -lextrudeModel \
-    -lcreateShellMesh \
     -lfiniteVolume \
     -lmeshTools \
-    -ldynamicMesh
+    -ldynamicMesh \
+    -lextrudeModel

applications/utilities/mesh/generation/extrude2DMesh/Make/options

@@ -3,7 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -Iextrude2DMesh/lnInclude \
-    -I../extrude/extrudeModel/lnInclude
+    -I$(LIB_SRC)/mesh/extrudeModel/lnInclude
 
 EXE_LIBS = \
     -lsurfMesh \

applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options

@@ -4,7 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(FOAM_APP)/utilities/mesh/generation/extrude/extrudeModel/lnInclude
+    -I$(LIB_SRC)/mesh/extrudeModel/lnInclude
 
 LIB_LIBS = \
     -lmeshTools \

applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C

@@ -257,7 +257,7 @@ int main(int argc, char *argv[])
         const pointField& points = mesh().points();
 
         const boundBox& bb = mesh().bounds();
-        const scalar mergeDim = 1E-4 * bb.minDim();
+        const scalar mergeDim = 1e-4 * bb.minDim();
 
         forAll(edges, edgeI)
         {

applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C

@@ -41,8 +41,8 @@ Foam::label Foam::findOppositeWedge
             if
             (
                 pp.size() == wpp.size()
-             && mag(pp.axis() & wpp.axis()) >= (1-1E-3)
-             && mag(ppCosAngle - wppCosAngle) >= 1E-3
+             && mag(pp.axis() & wpp.axis()) >= (1-1e-3)
+             && mag(ppCosAngle - wppCosAngle) >= 1e-3
             )
             {
                 return patchI;
@@ -106,7 +106,7 @@ bool Foam::checkWedges
             );
 
 
-            if (mag(opp.axis() & pp.axis()) < (1-1E-3))
+            if (mag(opp.axis() & pp.axis()) < (1-1e-3))
             {
                 if (report)
                 {

applications/utilities/mesh/manipulation/mirrorMesh/mirrorMesh.C

@@ -37,15 +37,24 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
+
+#   include "addOverwriteOption.H"
 #   include "setRootCase.H"
 #   include "createTime.H"
 
+    const bool overwrite = args.optionFound("overwrite");
+
+    if (!overwrite)
+    {
+        runTime++;
+    }
+
     mirrorFvMesh mesh
     (
         IOobject
         (
             mirrorFvMesh::defaultRegion,
-            runTime.timeName(),
+            runTime.constant(),
             runTime
         )
     );

applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.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
@@ -52,7 +52,7 @@ void Foam::meshDualiser::checkPolyTopoChange(const polyTopoChange& meshMod)
     label nUnique = mergePoints
     (
         points,
-        1E-6,
+        1e-6,
         false,
         oldToNew
     );
@@ -226,7 +226,7 @@ Foam::label Foam::meshDualiser::addInternalFace
         label nUnique = mergePoints
         (
             facePoints,
-            1E-6,
+            1e-6,
             false,
             oldToNew
         );

applications/utilities/mesh/manipulation/refineMesh/refineMesh.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
@@ -52,7 +52,7 @@ using namespace Foam;
 
 
 // Max cos angle for edges to be considered aligned with axis.
-static const scalar edgeTol = 1E-3;
+static const scalar edgeTol = 1e-3;
 
 
 // Calculate some edge statistics on mesh.
@@ -206,7 +206,7 @@ label twoDNess(const polyMesh& mesh)
             minLen = min(minLen, mesh.edges()[cEdges[i]].mag(mesh.points()));
         }
 
-        if (cellPlane.distance(ctrs[cellI]) > 1E-6*minLen)
+        if (cellPlane.distance(ctrs[cellI]) > 1e-6*minLen)
         {
             // Centres not in plane
             return  -1;
@@ -274,7 +274,7 @@ label twoDNess(const polyMesh& mesh)
 
             const scalarField cosAngle(mag(n/mag(n) & cellPlane.normal()));
 
-            if (mag(min(cosAngle) - max(cosAngle)) > 1E-6)
+            if (mag(min(cosAngle) - max(cosAngle)) > 1e-6)
             {
                 // cosAngle should be either ~1 over all faces (2D front and
                 // back) or ~0 (all other patches perp to 2D)

applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C

@@ -1036,6 +1036,22 @@ int main(int argc, char *argv[])
         (
             UIndirectList<label>(faceProcAddressing, map().faceMap())
         );
+
+        // Detect any flips.
+        const labelHashSet& fff = map().flipFaceFlux();
+        forAllConstIter(labelHashSet, fff, iter)
+        {
+            label faceI = iter.key();
+            label masterFaceI = faceProcAddressing[faceI];
+
+            faceProcAddressing[faceI] = -masterFaceI;
+
+            if (masterFaceI == 0)
+            {
+                FatalErrorIn(args.executable()) << "problem faceI:" << faceI
+                    << " masterFaceI:" << masterFaceI << exit(FatalError);
+            }
+        }
     }
     if (pointProcAddressing.headerOk())
     {
@@ -1083,9 +1099,13 @@ int main(int argc, char *argv[])
         );
         Info<< "After renumbering :" << nl
             << "    band           : " << band << nl
-            << "    profile        : " << profile << nl
-            << "    rms frontwidth : " << rmsFrontwidth << nl
-            << endl;
+            << "    profile        : " << profile << nl;
+        if (doFrontWidth)
+        {
+
+            Info<< "    rms frontwidth : " << rmsFrontwidth << nl;
+        }
+        Info<< endl;
     }
 
     if (orderPoints)

applications/utilities/parallelProcessing/decomposePar/decomposePar.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
@@ -43,7 +43,10 @@ Usage
     Copy any \a uniform directories too.
 
     \param -constant \n
-    Override controlDict settings and use constant directory.
+    \param -time xxx:yyy \n
+    Override controlDict settings and decompose selected times. Does not
+    re-decompose the mesh i.e. does not handle moving mesh or changing
+    mesh cases.
 
     \param -fields \n
     Use existing geometry decomposition and convert fields only.
@@ -123,11 +126,9 @@ int main(int argc, char *argv[])
         "ifRequired",
         "only decompose geometry if the number of domains has changed"
     );
-    argList::addBoolOption
-    (
-        "constant",
-        "include the 'constant/' dir in the times list"
-    );
+
+    // Include explicit constant options, have zero from time range
+    timeSelector::addOptions(true, false);
 
     #include "setRootCase.H"
 
@@ -146,23 +147,10 @@ int main(int argc, char *argv[])
     bool forceOverwrite          = args.optionFound("force");
     bool ifRequiredDecomposition = args.optionFound("ifRequired");
 
+    // Set time from database
     #include "createTime.H"
-
-    // Allow -constant to override controlDict settings.
-    if (args.optionFound("constant"))
-    {
-        instantList timeDirs = timeSelector::select0(runTime, args);
-        if (runTime.timeName() != runTime.constant())
-        {
-            FatalErrorIn(args.executable())
-                << "No '" << runTime.constant() << "' time present." << endl
-                << "Valid times are " << runTime.times()
-                << exit(FatalError);
-        }
-    }
-
-
-    Info<< "Time = " << runTime.timeName() << endl;
+    // Allow override of time
+    instantList times = timeSelector::selectIfPresent(runTime, args);
 
     // determine the existing processor count directly
     label nProcs = 0;
@@ -334,23 +322,44 @@ int main(int argc, char *argv[])
     }
 
 
+
+    // Caches
+    // ~~~~~~
+    // Cached processor meshes and maps. These are only preserved if running
+    // with multiple times.
+    PtrList<Time> processorDbList(mesh.nProcs());
+    PtrList<fvMesh> procMeshList(mesh.nProcs());
+    PtrList<labelIOList> faceProcAddressingList(mesh.nProcs());
+    PtrList<labelIOList> cellProcAddressingList(mesh.nProcs());
+    PtrList<labelIOList> boundaryProcAddressingList(mesh.nProcs());
+    PtrList<fvFieldDecomposer> fieldDecomposerList(mesh.nProcs());
+    PtrList<dimFieldDecomposer> dimFieldDecomposerList(mesh.nProcs());
+    PtrList<labelIOList> pointProcAddressingList(mesh.nProcs());
+    PtrList<pointFieldDecomposer> pointFieldDecomposerList(mesh.nProcs());
+
+
+
+    // Loop over all times
+    forAll(times, timeI)
+    {
+        runTime.setTime(times[timeI], timeI);
+
+        Info<< "Time = " << runTime.timeName() << endl;
+
         // Search for list of objects for this time
         IOobjectList objects(mesh, runTime.timeName());
 
+
         // Construct the vol fields
         // ~~~~~~~~~~~~~~~~~~~~~~~~
         PtrList<volScalarField> volScalarFields;
         readFields(mesh, objects, volScalarFields);
-
         PtrList<volVectorField> volVectorFields;
         readFields(mesh, objects, volVectorFields);
-
         PtrList<volSphericalTensorField> volSphericalTensorFields;
         readFields(mesh, objects, volSphericalTensorFields);
-
         PtrList<volSymmTensorField> volSymmTensorFields;
         readFields(mesh, objects, volSymmTensorFields);
-
         PtrList<volTensorField> volTensorFields;
         readFields(mesh, objects, volTensorFields);
 
@@ -359,17 +368,13 @@ int main(int argc, char *argv[])
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         PtrList<DimensionedField<scalar, volMesh> > dimScalarFields;
         readFields(mesh, objects, dimScalarFields);
-
         PtrList<DimensionedField<vector, volMesh> > dimVectorFields;
         readFields(mesh, objects, dimVectorFields);
-
         PtrList<DimensionedField<sphericalTensor, volMesh> >
             dimSphericalTensorFields;
         readFields(mesh, objects, dimSphericalTensorFields);
-
         PtrList<DimensionedField<symmTensor, volMesh> > dimSymmTensorFields;
         readFields(mesh, objects, dimSymmTensorFields);
-
         PtrList<DimensionedField<tensor, volMesh> > dimTensorFields;
         readFields(mesh, objects, dimTensorFields);
 
@@ -394,16 +399,12 @@ int main(int argc, char *argv[])
 
         PtrList<pointScalarField> pointScalarFields;
         readFields(pMesh, objects, pointScalarFields);
-
         PtrList<pointVectorField> pointVectorFields;
         readFields(pMesh, objects, pointVectorFields);
-
         PtrList<pointSphericalTensorField> pointSphericalTensorFields;
         readFields(pMesh, objects, pointSphericalTensorFields);
-
         PtrList<pointSymmTensorField> pointSymmTensorFields;
         readFields(pMesh, objects, pointSymmTensorFields);
-
         PtrList<pointTensorField> pointTensorFields;
         readFields(pMesh, objects, pointTensorFields);
 
@@ -419,24 +420,37 @@ int main(int argc, char *argv[])
         // Particles
         PtrList<Cloud<indexedParticle> > lagrangianPositions(cloudDirs.size());
         // Particles per cell
-    PtrList< List<SLList<indexedParticle*>*> > cellParticles(cloudDirs.size());
+        PtrList< List<SLList<indexedParticle*>*> > cellParticles
+        (
+            cloudDirs.size()
+        );
 
-    PtrList<PtrList<labelIOField> > lagrangianLabelFields(cloudDirs.size());
+        PtrList<PtrList<labelIOField> > lagrangianLabelFields
+        (
+            cloudDirs.size()
+        );
         PtrList<PtrList<labelFieldCompactIOField> > lagrangianLabelFieldFields
         (
             cloudDirs.size()
         );
-    PtrList<PtrList<scalarIOField> > lagrangianScalarFields(cloudDirs.size());
+        PtrList<PtrList<scalarIOField> > lagrangianScalarFields
+        (
+            cloudDirs.size()
+        );
         PtrList<PtrList<scalarFieldCompactIOField> > lagrangianScalarFieldFields
         (
             cloudDirs.size()
         );
-    PtrList<PtrList<vectorIOField> > lagrangianVectorFields(cloudDirs.size());
+        PtrList<PtrList<vectorIOField> > lagrangianVectorFields
+        (
+            cloudDirs.size()
+        );
         PtrList<PtrList<vectorFieldCompactIOField> > lagrangianVectorFieldFields
         (
             cloudDirs.size()
         );
-    PtrList<PtrList<sphericalTensorIOField> > lagrangianSphericalTensorFields
+        PtrList<PtrList<sphericalTensorIOField> >
+        lagrangianSphericalTensorFields
         (
             cloudDirs.size()
         );
@@ -478,7 +492,8 @@ int main(int argc, char *argv[])
                 // Read lagrangian particles
                 // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
-            Info<< "Identified lagrangian data set: " << cloudDirs[i] << endl;
+                Info<< "Identified lagrangian data set: " << cloudDirs[i]
+                    << endl;
 
                 lagrangianPositions.set
                 (
@@ -534,8 +549,8 @@ int main(int argc, char *argv[])
 
                     if (!cellParticles[cloudI][celli])
                     {
-                    cellParticles[cloudI][celli] = new SLList<indexedParticle*>
-                    ();
+                        cellParticles[cloudI][celli] =
+                            new SLList<indexedParticle*>();
                     }
 
                     cellParticles[cloudI][celli]->append(&iter());
@@ -676,13 +691,24 @@ int main(int argc, char *argv[])
         {
             Info<< "Processor " << procI << ": field transfer" << endl;
 
+
             // open the database
-        Time processorDb
+            if (!processorDbList.set(procI))
+            {
+                processorDbList.set
+                (
+                    procI,
+                    new Time
                     (
                         Time::controlDictName,
                         args.rootPath(),
-            args.caseName()/fileName(word("processor") + name(procI))
+                        args.caseName()
+                       /fileName(word("processor") + name(procI))
+                    )
                 );
+            }
+            Time& processorDb = processorDbList[procI];
+
 
             processorDb.setTime(runTime);
 
@@ -697,7 +723,12 @@ int main(int argc, char *argv[])
             }
 
             // read the mesh
-        fvMesh procMesh
+            if (!procMeshList.set(procI))
+            {
+                procMeshList.set
+                (
+                    procI,
+                    new fvMesh
                     (
                         IOobject
                         (
@@ -705,9 +736,18 @@ int main(int argc, char *argv[])
                             processorDb.timeName(),
                             processorDb
                         )
+                    )
                 );
+            }
+            const fvMesh& procMesh = procMeshList[procI];
 
-        labelIOList faceProcAddressing
+
+            if (!faceProcAddressingList.set(procI))
+            {
+                faceProcAddressingList.set
+                (
+                    procI,
+                    new labelIOList
                     (
                         IOobject
                         (
@@ -718,9 +758,19 @@ int main(int argc, char *argv[])
                             IOobject::MUST_READ,
                             IOobject::NO_WRITE
                         )
+                    )
                 );
+            }
+            const labelIOList& faceProcAddressing =
+                faceProcAddressingList[procI];
 
-        labelIOList cellProcAddressing
+
+            if (!cellProcAddressingList.set(procI))
+            {
+                cellProcAddressingList.set
+                (
+                    procI,
+                    new labelIOList
                     (
                         IOobject
                         (
@@ -731,9 +781,19 @@ int main(int argc, char *argv[])
                             IOobject::MUST_READ,
                             IOobject::NO_WRITE
                         )
+                    )
                 );
+            }
+            const labelIOList& cellProcAddressing =
+                cellProcAddressingList[procI];
 
-        labelIOList boundaryProcAddressing
+
+            if (!boundaryProcAddressingList.set(procI))
+            {
+                boundaryProcAddressingList.set
+                (
+                    procI,
+                    new labelIOList
                     (
                         IOobject
                         (
@@ -744,18 +804,32 @@ int main(int argc, char *argv[])
                             IOobject::MUST_READ,
                             IOobject::NO_WRITE
                         )
+                    )
                 );
+            }
+            const labelIOList& boundaryProcAddressing =
+                boundaryProcAddressingList[procI];
+
 
             // FV fields
             {
-            fvFieldDecomposer fieldDecomposer
+                if (!fieldDecomposerList.set(procI))
+                {
+                    fieldDecomposerList.set
+                    (
+                        procI,
+                        new fvFieldDecomposer
                         (
                             mesh,
                             procMesh,
                             faceProcAddressing,
                             cellProcAddressing,
                             boundaryProcAddressing
+                        )
                     );
+                }
+                const fvFieldDecomposer& fieldDecomposer =
+                    fieldDecomposerList[procI];
 
                 fieldDecomposer.decomposeFields(volScalarFields);
                 fieldDecomposer.decomposeFields(volVectorFields);
@@ -768,23 +842,43 @@ int main(int argc, char *argv[])
                 fieldDecomposer.decomposeFields(surfaceSphericalTensorFields);
                 fieldDecomposer.decomposeFields(surfaceSymmTensorFields);
                 fieldDecomposer.decomposeFields(surfaceTensorFields);
+
+                if (times.size() == 1)
+                {
+                    // Clear cached decomposer
+                    fieldDecomposerList.set(procI, NULL);
+                }
             }
 
             // Dimensioned fields
             {
-            dimFieldDecomposer fieldDecomposer
+                if (!dimFieldDecomposerList.set(procI))
+                {
+                    dimFieldDecomposerList.set
+                    (
+                        procI,
+                        new dimFieldDecomposer
                         (
                             mesh,
                             procMesh,
                             faceProcAddressing,
                             cellProcAddressing
+                        )
                     );
+                }
+                const dimFieldDecomposer& dimDecomposer =
+                    dimFieldDecomposerList[procI];
 
-            fieldDecomposer.decomposeFields(dimScalarFields);
-            fieldDecomposer.decomposeFields(dimVectorFields);
-            fieldDecomposer.decomposeFields(dimSphericalTensorFields);
-            fieldDecomposer.decomposeFields(dimSymmTensorFields);
-            fieldDecomposer.decomposeFields(dimTensorFields);
+                dimDecomposer.decomposeFields(dimScalarFields);
+                dimDecomposer.decomposeFields(dimVectorFields);
+                dimDecomposer.decomposeFields(dimSphericalTensorFields);
+                dimDecomposer.decomposeFields(dimSymmTensorFields);
+                dimDecomposer.decomposeFields(dimTensorFields);
+
+                if (times.size() == 1)
+                {
+                    dimFieldDecomposerList.set(procI, NULL);
+                }
             }
 
 
@@ -798,7 +892,12 @@ int main(int argc, char *argv[])
              || pointTensorFields.size()
             )
             {
-            labelIOList pointProcAddressing
+                if (!pointProcAddressingList.set(procI))
+                {
+                    pointProcAddressingList.set
+                    (
+                        procI,
+                        new labelIOList
                         (
                             IOobject
                             (
@@ -809,23 +908,43 @@ int main(int argc, char *argv[])
                                 IOobject::MUST_READ,
                                 IOobject::NO_WRITE
                             )
+                        )
                     );
+                }
+                const labelIOList& pointProcAddressing =
+                    pointProcAddressingList[procI];
 
-            pointMesh procPMesh(procMesh);
+                const pointMesh& procPMesh = pointMesh::New(procMesh);
 
-            pointFieldDecomposer fieldDecomposer
+                if (!pointFieldDecomposerList.set(procI))
+                {
+                    pointFieldDecomposerList.set
+                    (
+                        procI,
+                        new pointFieldDecomposer
                         (
                             pMesh,
                             procPMesh,
                             pointProcAddressing,
                             boundaryProcAddressing
+                        )
                     );
+                }
+                const pointFieldDecomposer& pointDecomposer =
+                    pointFieldDecomposerList[procI];
 
-            fieldDecomposer.decomposeFields(pointScalarFields);
-            fieldDecomposer.decomposeFields(pointVectorFields);
-            fieldDecomposer.decomposeFields(pointSphericalTensorFields);
-            fieldDecomposer.decomposeFields(pointSymmTensorFields);
-            fieldDecomposer.decomposeFields(pointTensorFields);
+                pointDecomposer.decomposeFields(pointScalarFields);
+                pointDecomposer.decomposeFields(pointVectorFields);
+                pointDecomposer.decomposeFields(pointSphericalTensorFields);
+                pointDecomposer.decomposeFields(pointSymmTensorFields);
+                pointDecomposer.decomposeFields(pointTensorFields);
+
+
+                if (times.size() == 1)
+                {
+                    pointProcAddressingList.set(procI, NULL);
+                    pointFieldDecomposerList.set(procI, NULL);
+                }
             }
 
 
@@ -940,9 +1059,23 @@ int main(int argc, char *argv[])
                     chDir(currentDir);
                 }
             }
-    }
 
 
+
+            // We have cached all the constant mesh data for the current
+            // processor. This is only important if running with multiple
+            // times, otherwise it is just extra storage.
+            if (times.size() == 1)
+            {
+                boundaryProcAddressingList.set(procI, NULL);
+                cellProcAddressingList.set(procI, NULL);
+                faceProcAddressingList.set(procI, NULL);
+                procMeshList.set(procI, NULL);
+                processorDbList.set(procI, NULL);
+            }
+        }
+    }
+
     Info<< "\nEnd.\n" << endl;
 
     return 0;

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.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
@@ -55,7 +55,7 @@ using namespace Foam;
 
 // Tolerance (as fraction of the bounding box). Needs to be fairly lax since
 // usually meshes get written with limited precision (6 digits)
-static const scalar defaultMergeTol = 1E-7;
+static const scalar defaultMergeTol = 1e-7;
 
 
 static void renumber
@@ -290,7 +290,7 @@ int main(int argc, char *argv[])
         "mergeTol",
         "scalar",
         "specify the merge distance relative to the bounding box size "
-        "(default 1E-7)"
+        "(default 1e-7)"
     );
     argList::addBoolOption
     (

applications/utilities/parallelProcessing/redistributePar/redistributePar.C

@@ -60,7 +60,7 @@ Description
 
 // Tolerance (as fraction of the bounding box). Needs to be fairly lax since
 // usually meshes get written with limited precision (6 digits)
-static const scalar defaultMergeTol = 1E-6;
+static const scalar defaultMergeTol = 1e-6;
 
 
 //// Read mesh if available. Otherwise create empty mesh with same non-proc
@@ -721,7 +721,7 @@ int main(int argc, char *argv[])
         "mergeTol",
         "scalar",
         "specify the merge distance relative to the bounding box size "
-        "(default 1E-6)"
+        "(default 1e-6)"
     );
 #   include "setRootCase.H"
 

applications/utilities/postProcessing/miscellaneous/temporalInterpolate/temporalInterpolate.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
@@ -37,6 +37,8 @@ Description
 #include "surfaceFields.H"
 #include "pointFields.H"
 #include "ReadFields.H"
+#include "interpolationWeights.H"
+#include "uniformInterpolate.H"
 
 using namespace Foam;
 
@@ -48,6 +50,8 @@ class fieldInterpolator
     const HashSet<word>& selectedFields_;
     instant ti_;
     instant ti1_;
+    const interpolationWeights& interpolator_;
+    const wordList& timeNames_;
     int divisions_;
 
 public:
@@ -60,6 +64,8 @@ public:
         const HashSet<word>& selectedFields,
         const instant& ti,
         const instant& ti1,
+        const interpolationWeights& interpolator,
+        const wordList& timeNames,
         int divisions
     )
     :
@@ -69,6 +75,8 @@ public:
         selectedFields_(selectedFields),
         ti_(ti),
         ti1_(ti1),
+        interpolator_(interpolator),
+        timeNames_(timeNames),
         divisions_(divisions)
     {}
 
@@ -98,34 +106,6 @@ void fieldInterpolator::interpolate()
             {
                 Info<< " " << fieldIter()->name() << '(';
 
-                GeoFieldType fieldi
-                (
-                    IOobject
-                    (
-                        fieldIter()->name(),
-                        ti_.name(),
-                        fieldIter()->db(),
-                        IOobject::MUST_READ,
-                        IOobject::NO_WRITE,
-                        false
-                    ),
-                    mesh_
-                );
-
-                GeoFieldType fieldi1
-                (
-                    IOobject
-                    (
-                        fieldIter()->name(),
-                        ti1_.name(),
-                        fieldIter()->db(),
-                        IOobject::MUST_READ,
-                        IOobject::NO_WRITE,
-                        false
-                    ),
-                    mesh_
-                );
-
                 scalar deltaT = (ti1_.value() - ti_.value())/(divisions_ + 1);
 
                 for (int j=0; j<divisions_; j++)
@@ -141,20 +121,51 @@ void fieldInterpolator::interpolate()
                         Info<< " ";
                     }
 
-                    scalar lambda = scalar(j + 1)/scalar(divisions_ + 1);
+                    // Calculate times to read and weights
+                    labelList indices;
+                    scalarField weights;
+                    interpolator_.valueWeights
+                    (
+                        runTime_.value(),
+                        indices,
+                        weights
+                    );
+
+                    const wordList selectedTimeNames
+                    (
+                        UIndirectList<word>(timeNames_, indices)()
+                    );
+
+                    //Info<< "For time " << runTime_.value()
+                    //    << " need times " << selectedTimeNames
+                    //    << " need weights " << weights << endl;
+
+
+                    // Read on the objectRegistry all the required fields
+                    ReadFields<GeoFieldType>
+                    (
+                        fieldIter()->name(),
+                        mesh_,
+                        selectedTimeNames
+                    );
 
                     GeoFieldType fieldj
+                    (
+                        uniformInterpolate<GeoFieldType>
                         (
                             IOobject
                             (
                                 fieldIter()->name(),
-                            timej.name(),
+                                runTime_.timeName(),
                                 fieldIter()->db(),
                                 IOobject::NO_READ,
                                 IOobject::NO_WRITE,
                                 false
                             ),
-                        (1.0 - lambda)*fieldi + lambda*fieldi1
+                            fieldIter()->name(),
+                            selectedTimeNames,
+                            weights
+                        )
                     );
 
                     fieldj.write();
@@ -188,6 +199,12 @@ int main(int argc, char *argv[])
         "integer",
         "specify number of temporal sub-divisions to create (default = 1)."
     );
+    argList::addOption
+    (
+        "interpolationType",
+        "word",
+        "specify type of interpolation (linear or spline)"
+    );
 
     #include "setRootCase.H"
     #include "createTime.H"
@@ -198,15 +215,51 @@ int main(int argc, char *argv[])
     {
         args.optionLookup("fields")() >> selectedFields;
     }
+    if (selectedFields.size())
+    {
+        Info<< "Interpolating fields " << selectedFields << nl << endl;
+    }
+    else
+    {
+        Info<< "Interpolating all fields" << nl << endl;
+    }
+
 
     int divisions = 1;
     if (args.optionFound("divisions"))
     {
         args.optionLookup("divisions")() >> divisions;
     }
+    Info<< "Using " << divisions << " per time interval" << nl << endl;
+
+
+    const word interpolationType = args.optionLookupOrDefault<word>
+    (
+        "interpolationType",
+        "linear"
+    );
+    Info<< "Using interpolation " << interpolationType << nl << endl;
+
 
     instantList timeDirs = timeSelector::select0(runTime, args);
 
+    scalarField timeVals(timeDirs.size());
+    wordList timeNames(timeDirs.size());
+    forAll(timeDirs, i)
+    {
+        timeVals[i] = timeDirs[i].value();
+        timeNames[i] = timeDirs[i].name();
+    }
+    autoPtr<interpolationWeights> interpolatorPtr
+    (
+        interpolationWeights::New
+        (
+            interpolationType,
+            timeVals
+        )
+    );
+
+
     #include "createMesh.H"
 
     Info<< "Interpolating fields for times:" << endl;
@@ -226,6 +279,8 @@ int main(int argc, char *argv[])
             selectedFields,
             timeDirs[timei],
             timeDirs[timei+1],
+            interpolatorPtr(),
+            timeNames,
             divisions
         );
 

applications/utilities/preProcessing/changeDictionary/changeDictionary.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
@@ -70,6 +70,7 @@ Usage
 #include "IOPtrList.H"
 #include "volFields.H"
 #include "stringListOps.H"
+#include "timeSelector.H"
 
 using namespace Foam;
 
@@ -256,8 +257,12 @@ int main(int argc, char *argv[])
     (
         "instance",
         "name",
-        "specify alternate time instance - default is latest time"
+        "override instance setting (default is the time name)"
     );
+
+    // Add explicit time option
+    timeSelector::addOptions();
+
     argList::addBoolOption
     (
         "literalRE",
@@ -272,6 +277,18 @@ int main(int argc, char *argv[])
 
     #include "setRootCase.H"
     #include "createTime.H"
+
+    // Optionally override controlDict time with -time options
+    instantList times = timeSelector::selectIfPresent(runTime, args);
+    if (times.size() < 1)
+    {
+        FatalErrorIn(args.executable())
+            << "No times selected." << exit(FatalError);
+    }
+    runTime.setTime(times[0], 0);
+    word instance = args.optionLookupOrDefault("instance", runTime.timeName());
+
+
     #include "createNamedMesh.H"
 
     const word dictName("changeDictionaryDict");
@@ -317,11 +334,6 @@ int main(int argc, char *argv[])
         regionPrefix = regionName;
     }
 
-    word instance = runTime.timeName();
-    if (args.options().found("instance"))
-    {
-        instance = args.options()["instance"];
-    }
 
     // Make sure we do not use the master-only reading since we read
     // fields (different per processor) as dictionaries.

applications/utilities/preProcessing/mapFields/mapLagrangian.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
@@ -32,7 +32,7 @@ License
 namespace Foam
 {
 
-static const scalar perturbFactor = 1E-6;
+static const scalar perturbFactor = 1e-6;
 
 
 // Special version of findCell that generates a cell guaranteed to be

applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H

@@ -7,7 +7,7 @@ List<treeBoundBox> meshBb
     treeBoundBox
     (
         boundBox(coarseMesh.points(), false)
-    ).extend(rndGen, 1E-3)
+    ).extend(rndGen, 1e-3)
 );
 
 // Dummy bounds dictionary

applications/utilities/surface/surfaceCheck/surfaceCheck.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
@@ -367,7 +367,7 @@ int main(int argc, char *argv[])
 
             forAll(triQ, faceI)
             {
-                if (triQ[faceI] < 1E-11)
+                if (triQ[faceI] < 1e-11)
                 {
                     problemFaces.append(faceI);
                 }
@@ -427,9 +427,9 @@ int main(int argc, char *argv[])
         const pointField& localPoints = surf.localPoints();
 
         const boundBox bb(localPoints);
-        scalar smallDim = 1E-6 * bb.mag();
+        scalar smallDim = 1e-6 * bb.mag();
 
-        Info<< "Checking for points less than 1E-6 of bounding box ("
+        Info<< "Checking for points less than 1e-6 of bounding box ("
             << bb.span() << " meter) apart."
             << endl;
 

applications/utilities/surface/surfaceFeatureExtract/Allwmake

@@ -15,7 +15,9 @@ then
            CGAL_LIBDIR=-L$CGAL_ARCH_PATH/lib \
            LAPACK_LIB=-llapack \
            BLAS_LIB=-lblas \
-           CGAL_LIB=-lCGAL"
+           CGAL_LIB=-lCGAL \
+           CGAL_BOOST_LIB=-lboost_thread \
+           CGAL_MPFR_LIB=-lmpfr"
 else
     echo
     echo "Compiling surfaceFeatureExtract without CGAL curvature support"

applications/utilities/surface/surfaceFeatureExtract/Make/options

@@ -15,6 +15,8 @@ EXE_INC = \
 
 EXE_LIBS = \
     $(CGAL_LIBS) \
+    ${CGAL_BOOST_LIB} \
+    ${CGAL_MPFR_LIB} \
     ${CGAL_LIBDIR} \
     ${LAPACK_LIB} \
     ${BLAS_LIB} \

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict

@@ -30,9 +30,6 @@ surface1.stl
 
     // Write options
 
-        // Write .eMesh file (for snappyHexMesh)
-        writeFeatureEdgeMesh    yes;
-
         // Write features to obj format for postprocessing
         writeObj                yes;
 }

applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.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
@@ -159,7 +159,7 @@ int main(int argc, char *argv[])
             treeBoundBox
             (
                 boundBox(mesh.points(), false)
-            ).extend(rndGen, 1E-3)
+            ).extend(rndGen, 1e-3)
         );
         Pstream::gatherList(meshBb);
         Pstream::scatterList(meshBb);

applications/utilities/surface/surfaceSubset/surfaceSubset.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
@@ -254,7 +254,7 @@ int main(int argc, char *argv[])
                 selectSurf,
                 indexedOctree<treeDataTriSurface>::perturbTol()
             ),
-            bb.extend(rndGen, 1E-4),    // slightly randomize bb
+            bb.extend(rndGen, 1e-4),    // slightly randomize bb
             8,      // maxLevel
             10,     // leafsize
             3.0     // duplicity