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

applications/solvers/combustion/fireFoam/pEqn.H

@@ -5,7 +5,7 @@ surfaceScalarField rhorAUf(rAU.name() + 'f', 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

@@ -5,14 +5,14 @@
     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

@@ -6,7 +6,7 @@
     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/buoyantSimpleFoam/pEqn.H

@@ -9,7 +9,7 @@
     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/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -5,7 +5,7 @@ surfaceScalarField rhorAUf(rAU.name() + 'f', 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/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/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 field T\n" << endl;
-    Info<< "Reading transportProperties\n" << endl;
+    volScalarField T
-    twoPhaseMixture twoPhaseProperties(U, phi);
-
-    dimensionedScalar rho10
     (
-        twoPhaseProperties.subDict
+        IOobject
         (
-            twoPhaseProperties.phase1Name()
+            "T",
-        ).lookup("rho0")
+            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
+    (
+        "alpha" + twoPhaseProperties.phase2Name(),
+        scalar(1) - alpha1
+    );
+
+    dimensionedScalar k1
+    (
+        "k",
+        dimensionSet(1, 1, -3, -1, 0),
+        twoPhaseProperties.subDict
         (
-            (p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
+            twoPhaseProperties.phase1Name()
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
+        ).lookup("k")
-            pMin
+    );
+
+    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);
+    autoPtr<phaseEquationOfState> eos2
-    volScalarField rho2(rho20 + psi2*p);
+    (
+        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());
+    rho1 = eos1->rho(p, T);
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    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 =
+        surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
-        (
+        surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
-            fvm::ddt(p_rgh)
+
-          + fvm::div(phi, p_rgh)
+        p_rghEqnComp1 =
-          - fvm::Sp(fvc::div(phi), p_rgh)
+            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
-        );
+          + fvc::div(phid1, p_rgh)
-    }
+          - fvc::Sp(fvc::div(phid1), p_rgh);
-    else
+
-    {
+        p_rghEqnComp2 =
-        p_rghEqnComp =
+            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
-        (
+          + fvc::div(phid2, p_rgh)
-            fvm::ddt(p_rgh)
+          - fvc::Sp(fvc::div(phid2), p_rgh);
-          + fvc::div(phi, p_rgh)
-          - fvc::Sp(fvc::div(phi), p_rgh)
-        );
     }
 
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    U = rAU*UEqn.H();
+    surfaceScalarField phiHbyA
-
-    surfaceScalarField phiU
     (
-        "phiU",
+        "phiHbyA",
-        (fvc::interpolate(U) & mesh.Sf())
+        (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(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
-              + max(alpha2, scalar(0))*(psi2/rho2)
+              + (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);
+                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
-            phi += p_rghEqnIncomp.flux();
+              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
+            );
+
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+    p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-    U.correctBoundaryConditions();
 
-    p = max
+    rho1 = eos1->rho(p, T);
-    (
+    rho2 = eos2->rho(p, T);
-        (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;

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());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
-    surfaceScalarField phiU
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-        "phiU",
+        "phiHbyA",
-        (fvc::interpolate(U) & mesh.Sf())
+        (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)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-    )*rAUf*mesh.magSf();
+          - 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));
@@ -32,13 +38,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;

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,38 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    phiAbs =
+    surfaceScalarField phiHbyA
-        (fvc::interpolate(U) & mesh.Sf())
+    (
-      + fvc::ddtPhiCorr(rAU, rho, U, phiAbs);
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+    );
 
     if (p_rgh.needReference())
     {
-        fvc::makeRelative(phiAbs, U);
+        fvc::makeRelative(phiHbyA, U);
-        adjustPhi(phiAbs, U, p_rgh);
+        adjustPhi(phiHbyA, U, p_rgh);
-        fvc::makeAbsolute(phiAbs, U);
+        fvc::makeAbsolute(phiHbyA, U);
     }
 
-    phi = phiAbs +
+    surfaceScalarField phig
     (
-        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+        (
-      - ghf*fvc::snGrad(rho)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-    )*rAUf*mesh.magSf();
+          - 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,13 +41,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiAbs)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     phiAbs = phi;

applications/solvers/multiphase/interFoam/pEqn.H

@@ -1,28 +1,34 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
-    surfaceScalarField phiU
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-        "phiU",
+        "phiHbyA",
-        (fvc::interpolate(U) & mesh.Sf())
+        (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)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-    )*rAUf*mesh.magSf();
+          - 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));
@@ -31,13 +37,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -1,24 +1,28 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
+        "phiHbyA",
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi =
+    surfaceScalarField phig
-        phiU
+    (
-      + (
+        (
             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,13 +42,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += p_rghEqn.flux();
+            phi = phiHbyA + p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;

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

@@ -1,31 +1,32 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
+        "phiHbyA",
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
-      //+ fvc::ddtPhiCorr(rAU, rho, U, phi)
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiU);
+    mrfZones.relativeFlux(phiHbyA);
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi =
+    surfaceScalarField phig
-        phiU
+    (
-      + (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
-            mixture.surfaceTensionForce()
+    );
-          - 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,13 +35,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;

applications/solvers/multiphase/multiphaseInterFoam/pEqn.H

@@ -1,29 +1,34 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
+        "phiHbyA",
-        (fvc::interpolate(U) & mesh.Sf())
+        (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)
+            mixture.surfaceTensionForce()
-    )*rAUf*mesh.magSf();
+          - 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));
@@ -32,13 +37,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;

applications/solvers/multiphase/settlingFoam/pEqn.H

@@ -1,54 +1,60 @@
-volScalarField rAU(1.0/UEqn.A());
-
-surfaceScalarField rAUf
-(
-    "(rho*(1|A(U)))",
-    fvc::interpolate(rho)*fvc::interpolate(rAU)
-);
-
-U = rAU*UEqn.H();
-phi =
-    fvc::interpolate(rho)
-   *(
-       (fvc::interpolate(U) & mesh.Sf())
-     + fvc::ddtPhiCorr(rAU, rho, U, phi)
-    );
-
-surfaceScalarField phiU("phiU", phi);
-phi -= ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
-
-while (pimple.correctNonOrthogonal())
 {
-    fvScalarMatrix p_rghEqn
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rho*rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-        fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phi)
+        fvc::interpolate(rho)
+       *(
+           (fvc::interpolate(HbyA) & mesh.Sf())
+         + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        )
     );
 
-    p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );
 
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+    phiHbyA += phig;
 
-    if (pimple.finalNonOrthogonalIter())
+    while (pimple.correctNonOrthogonal())
     {
-        phi -= p_rghEqn.flux();
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+        }
     }
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
 }
-
-p == p_rgh + rho*gh;
-
-if (p_rgh.needReference())
-{
-    p += dimensionedScalar
-    (
-        "p",
-        p.dimensions(),
-        pRefValue - getRefCellValue(p, pRefCell)
-    );
-    p_rgh = p - rho*gh;
-}
-
-#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());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
-    surfaceScalarField phiU
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-        "phiU",
+        "phiHbyA",
-        (fvc::interpolate(U) & mesh.Sf())
+        (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,13 +34,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;

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/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/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/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/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/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(pp.axis() & wpp.axis()) >= (1-1e-3)
-             && mag(ppCosAngle - wppCosAngle) >= 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/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/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
                     (
-                        IOobject
+                        uniformInterpolate<GeoFieldType>
                         (
+                            IOobject
+                            (
+                                fieldIter()->name(),
+                                runTime_.timeName(),
+                                fieldIter()->db(),
+                                IOobject::NO_READ,
+                                IOobject::NO_WRITE,
+                                false
+                            ),
                             fieldIter()->name(),
-                            timej.name(),
+                            selectedTimeNames,
-                            fieldIter()->db(),
+                            weights
-                            IOobject::NO_READ,
+                        )
-                            IOobject::NO_WRITE,
-                            false
-                        ),
-                        (1.0 - lambda)*fieldi + lambda*fieldi1
                     );
 
                     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;
 
@@ -252,12 +253,10 @@ int main(int argc, char *argv[])
         "file",
         "specify an alternative to system/changeDictionaryDict"
     );
-    argList::addOption
+
-    (
+    // Add explicit time option
-        "instance",
+    timeSelector::addOptions();
-        "name",
+
-        "specify alternate time instance - default is latest time"
-    );
     argList::addBoolOption
     (
         "literalRE",
@@ -272,6 +271,17 @@ 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);
+
+
     #include "createNamedMesh.H"
 
     const word dictName("changeDictionaryDict");
@@ -317,11 +327,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.
@@ -460,7 +465,7 @@ int main(int argc, char *argv[])
                 IOobject
                 (
                     fieldName,
-                    instance,
+                    runTime.timeName(),
                     mesh,
                     IOobject::MUST_READ_IF_MODIFIED,
                     IOobject::NO_WRITE,

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/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

src/OSspecific/POSIX/clockTime/clockTime.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
@@ -36,7 +36,7 @@ void Foam::clockTime::getTime(timeType& t)
 
 double Foam::clockTime::timeDifference(const timeType& beg, const timeType& end)
 {
-    return end.tv_sec - beg.tv_sec + 1E-6*(end.tv_usec - beg.tv_usec);
+    return end.tv_sec - beg.tv_sec + 1e-6*(end.tv_usec - beg.tv_usec);
 }
 
 

src/OpenFOAM/Make/files

@@ -576,6 +576,13 @@ $(interpolations)/interpolationTable/tableReaders/tableReaders.C
 $(interpolations)/interpolationTable/tableReaders/openFoam/openFoamTableReaders.C
 $(interpolations)/interpolationTable/tableReaders/csv/csvTableReaders.C
 
+interpolationWeights = $(interpolations)/interpolationWeights
+$(interpolationWeights)/interpolationWeights/interpolationWeights.C
+$(interpolationWeights)/linearInterpolationWeights/linearInterpolationWeights.C
+$(interpolationWeights)/splineInterpolationWeights/splineInterpolationWeights.C
+
+
+
 algorithms/indexedOctree/indexedOctreeName.C
 algorithms/indexedOctree/treeDataCell.C
 

src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H

@@ -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
@@ -119,6 +119,12 @@ public:
             inline void operator=(const T&);
 
 
+    // STL type definitions
+
+        //- Type of values the UList contains.
+        typedef T value_type;
+
+
     // Ostream operator
 
         //- Write UIndirectList to Ostream

src/OpenFOAM/db/objectRegistry/objectRegistry.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
@@ -151,9 +151,25 @@ Foam::wordList Foam::objectRegistry::sortedNames(const word& ClassName) const
 
 const Foam::objectRegistry& Foam::objectRegistry::subRegistry
 (
-    const word& name
+    const word& name,
+    const bool forceCreate
 ) const
 {
+    if (forceCreate && !foundObject<objectRegistry>(name))
+    {
+        objectRegistry* fieldsCachePtr = new objectRegistry
+        (
+            IOobject
+            (
+                name,
+                time().constant(),
+                *this,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            )
+        );
+        fieldsCachePtr->store();
+    }
     return lookupObject<objectRegistry>(name);
 }
 

src/OpenFOAM/db/objectRegistry/objectRegistry.H

@@ -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
@@ -147,8 +147,13 @@ public:
             template<class Type>
             wordList names() const;
 
-            //- Lookup and return a const sub-objectRegistry
+            //- Lookup and return a const sub-objectRegistry. Optionally create
-            const objectRegistry& subRegistry(const word& name) const;
+            //  it if it does not exist.
+            const objectRegistry& subRegistry
+            (
+                const word& name,
+                const bool forceCreate = false
+            ) const;
 
             //- Lookup and return all objects of the given Type
             template<class Type>

src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.C

@@ -131,6 +131,16 @@ dimensioned<Type>::dimensioned
     token nextToken(is);
     is.putBack(nextToken);
 
+    // Check if the original format is used in which the name is provided
+    // and reset the name to that read
+    if (nextToken.isWord())
+    {
+        is >> name_;
+        is >> nextToken;
+        is.putBack(nextToken);
+    }
+
+    // If the dimensions are provided compare with the argument
     if (nextToken == token::BEGIN_SQR)
     {
         dimensionSet dims(is);
@@ -391,8 +401,26 @@ dimensioned<Type> min
 template <class Type>
 Istream& operator>>(Istream& is, dimensioned<Type>& dt)
 {
-    // do a stream read op for a Type and a dimensions()et
+    token nextToken(is);
-    is >> dt.name_ >> dt.dimensions_ >> dt.value_;
+    is.putBack(nextToken);
+
+    // Check if the original format is used in which the name is provided
+    // and reset the name to that read
+    if (nextToken.isWord())
+    {
+        is >> dt.name_;
+        is >> nextToken;
+        is.putBack(nextToken);
+    }
+
+    // If the dimensions are provided reset the dimensions to those read
+    if (nextToken == token::BEGIN_SQR)
+    {
+        is >> dt.dimensions_;
+    }
+
+    // Read the value
+    is >> dt.value_;
 
     // Check state of Istream
     is.check("Istream& operator>>(Istream&, dimensioned<Type>&)");

src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class GeoField>
+Foam::tmp<GeoField> Foam::uniformInterpolate
+(
+    const HashPtrTable<GeoField, label, Hash<label> >& fields,
+    const labelList& indices,
+    const scalarField& weights
+)
+{
+    const GeoField& field0 = *(*fields.begin());
+
+    // Interpolate
+    tmp<GeoField> tfld
+    (
+        new GeoField
+        (
+            IOobject
+            (
+                "uniformInterpolate(" + field0.name() + ')',
+                field0.time().timeName(),
+                field0.db(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            weights[0]*(*fields[indices[0]])
+        )
+    );
+    GeoField& fld = tfld();
+
+    for (label i = 1; i < indices.size(); ++i)
+    {
+        fld += weights[i]*(*fields[indices[i]]);
+    }
+
+    return tfld;
+}
+
+
+template<class GeoField>
+Foam::tmp<GeoField> Foam::uniformInterpolate
+(
+    const IOobject& fieldIO,
+    const word& fieldName,
+    const wordList& times,
+    const scalarField& weights,
+    const objectRegistry& fieldsCache
+)
+{
+    // Look up the first field
+    const objectRegistry& time0Fields = fieldsCache.lookupObject
+    <
+        const objectRegistry
+    >
+    (
+        times[0]
+    );
+    const GeoField& field0 = time0Fields.lookupObject
+    <
+        const GeoField
+    >
+    (
+        fieldName
+    );
+
+
+    // Interpolate
+    tmp<GeoField> tfld(new GeoField(fieldIO, weights[0]*field0));
+    GeoField& fld = tfld();
+
+    for (label i = 1; i < times.size(); ++i)
+    {
+        const objectRegistry& timeIFields = fieldsCache.lookupObject
+        <
+            const objectRegistry
+        >
+        (
+            times[i]
+        );
+        const GeoField& fieldI = timeIFields.lookupObject
+        <
+            const GeoField
+        >
+        (
+            fieldName
+        );
+
+        fld += weights[i]*fieldI;
+    }
+
+    return tfld;
+}
+
+
+template<class GeoField>
+Foam::tmp<GeoField> Foam::uniformInterpolate
+(
+    const IOobject& fieldIO,
+    const word& fieldName,
+    const wordList& times,
+    const scalarField& weights,
+    const word& registryName
+)
+{
+    return uniformInterpolate<GeoField>
+    (
+        fieldIO,
+        fieldName,
+        times,
+        weights,
+        fieldIO.db().subRegistry(registryName, true)
+    );
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H

@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "GeometricField.H"
+#include "HashPtrTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Global functions  * * * * * * * * * * * * * //
+
+//- Interpolate selected fields (given by indices and corresponding weights)
+//  (boundary type becomes calculated). Fields stored per index. Field gets name
+//  "uniformInterpolate(" + fld.name() + ')'.
+template<class GeoField>
+tmp<GeoField> uniformInterpolate
+(
+    const HashPtrTable<GeoField, label, Hash<label> >& fields,
+    const labelList& indices,
+    const scalarField& weights
+);
+
+//- Interpolate fields. fieldsCache contains per timeName all loaded fields.
+//  Resulting field gets properties according to fieldIO
+template<class GeoField>
+tmp<GeoField> uniformInterpolate
+(
+    const IOobject& fieldIO,
+    const word& fieldName,
+    const wordList& times,
+    const scalarField& weights,
+    const objectRegistry& fieldsCache
+);
+
+//- Interpolate fields. fieldsCache contains per timeName all loaded fields.
+//  Resulting field gets properties according to fieldIO
+template<class GeoField>
+tmp<GeoField> uniformInterpolate
+(
+    const IOobject& fieldIO,
+    const word& fieldName,
+    const wordList& times,
+    const scalarField& weights,
+    const word& registryName = "fieldsCache"
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "uniformInterpolate.C"
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/ReadFields/ReadFields.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
@@ -127,4 +127,134 @@ Foam::wordList Foam::ReadFields
 }
 
 
+template<class GeoField>
+void Foam::ReadFields
+(
+    const word& fieldName,
+    const typename GeoField::Mesh& mesh,
+    const wordList& timeNames,
+    objectRegistry& fieldsCache
+)
+{
+    // Collect all times that are no longer used
+    {
+        HashSet<word> usedTimes(timeNames);
+
+        DynamicList<word> unusedTimes(fieldsCache.size());
+
+        forAllIter(objectRegistry, fieldsCache, timeIter)
+        {
+            const word& tm = timeIter.key();
+            if (!usedTimes.found(tm))
+            {
+                unusedTimes.append(tm);
+            }
+        }
+
+        //Info<< "Unloading times " << unusedTimes << endl;
+
+        forAll(unusedTimes, i)
+        {
+            objectRegistry& timeCache = const_cast<objectRegistry&>
+            (
+                fieldsCache.lookupObject<objectRegistry>(unusedTimes[i])
+            );
+            fieldsCache.checkOut(timeCache);
+        }
+    }
+
+
+    // Load any new fields
+    forAll(timeNames, i)
+    {
+        const word& tm = timeNames[i];
+
+        // Create if not found
+        if (!fieldsCache.found(tm))
+        {
+            //Info<< "Creating registry for time " << tm << endl;
+
+            // Create objectRegistry if not found
+            objectRegistry* timeCachePtr = new objectRegistry
+            (
+                IOobject
+                (
+                    tm,
+                    tm,
+                    fieldsCache,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                )
+            );
+            timeCachePtr->store();
+        }
+
+        // Obtain cache for current time
+        const objectRegistry& timeCache =
+            fieldsCache.lookupObject<objectRegistry>
+            (
+                tm
+            );
+
+        // Store field if not found
+        if (!timeCache.found(fieldName))
+        {
+            //Info<< "Loading field " << fieldName
+            //    << " for time " << tm << endl;
+
+            GeoField loadedFld
+            (
+                IOobject
+                (
+                    fieldName,
+                    tm,
+                    mesh.thisDb(),
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh
+            );
+
+            // Transfer to timeCache (new objectRegistry and store flag)
+            GeoField* fldPtr = new GeoField
+            (
+                IOobject
+                (
+                    fieldName,
+                    tm,
+                    timeCache,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                loadedFld
+            );
+            fldPtr->store();
+        }
+    }
+}
+
+
+template<class GeoField>
+void Foam::ReadFields
+(
+    const word& fieldName,
+    const typename GeoField::Mesh& mesh,
+    const wordList& timeNames,
+    const word& registryName
+)
+{
+    ReadFields<GeoField>
+    (
+        fieldName,
+        mesh,
+        timeNames,
+        const_cast<objectRegistry&>
+        (
+            mesh.thisDb().subRegistry(registryName, true)
+        )
+    );
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/fields/ReadFields/ReadFields.H

@@ -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,6 +55,28 @@ wordList ReadFields
     const bool syncPar = true
 );
 
+//- Helper routine to read GeometricFields. The fieldsCache is per time
+//  an objectRegistry of all stored fields
+template<class GeoField>
+static void ReadFields
+(
+    const word& fieldName,
+    const typename GeoField::Mesh& mesh,
+    const wordList& timeNames,
+    objectRegistry& fieldsCache
+);
+
+//- Helper routine to read GeometricFields. The fieldsCache is per time
+//  an objectRegistry of all stored fields
+template<class GeoField>
+static void ReadFields
+(
+    const word& fieldName,
+    const typename GeoField::Mesh& mesh,
+    const wordList& timeNames,
+    const word& registryName = "fieldsCache"
+);
+
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "interpolationWeights.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(interpolationWeights, 0);
+defineRunTimeSelectionTable(interpolationWeights, word);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+interpolationWeights::interpolationWeights
+(
+    const scalarField& samples
+)
+:
+    samples_(samples)
+{}
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+autoPtr<interpolationWeights> interpolationWeights::New
+(
+    const word& type,
+    const scalarField& samples
+)
+{
+    Info<< nl << "Selecting interpolationWeights "
+        << type << endl;
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(type);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "interpolationWeights::New(const word&, "
+            "const scalarField&)"
+        )   << "Unknown interpolationWeights type "
+            << type
+            << endl << endl
+            << "Valid interpolationWeights types are :" << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<interpolationWeights>(cstrIter()(samples));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+interpolationWeights::~interpolationWeights()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//objectRegistry& interpolationWeights::registry
+//(
+//    const objectRegistry& obr,
+//    const word& name
+//)
+//{
+//    if (!obr.foundObject<objectRegistry>(name))
+//    {
+//        objectRegistry* fieldsCachePtr = new objectRegistry
+//        (
+//            IOobject
+//            (
+//                name,
+//                obr.time().constant(),
+//                obr,
+//                IOobject::NO_READ,
+//                IOobject::NO_WRITE
+//            )
+//        );
+//        fieldsCachePtr->store();
+//    }
+//    return const_cast<objectRegistry&>(obr.subRegistry(name));
+//}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::interpolationWeights
+
+Description
+    Abstract base class for interpolating in 1D
+
+SourceFiles
+    interpolationWeights.C
+    interpolationWeightsTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interpolationWeights_H
+#define interpolationWeights_H
+
+#include "scalarField.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+class objectRegistry;
+
+/*---------------------------------------------------------------------------*\
+                       Class interpolationWeights Declaration
+\*---------------------------------------------------------------------------*/
+
+class interpolationWeights
+{
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        interpolationWeights(const interpolationWeights&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const interpolationWeights&);
+
+protected:
+
+        const scalarField& samples_;
+
+public:
+
+    //- Runtime type information
+    TypeName("interpolationWeights");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            interpolationWeights,
+            word,
+            (
+                const scalarField& samples
+            ),
+            (samples)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        interpolationWeights(const scalarField& samples);
+
+
+    // Selectors
+
+        //- Return a reference to the selected interpolationWeights
+        static autoPtr<interpolationWeights> New
+        (
+            const word& type,
+            const scalarField& samples
+        );
+
+
+    //- Destructor
+    virtual ~interpolationWeights();
+
+
+    // Member Functions
+
+        //- Calculate weights and indices to calculate t from samples.
+        //  Returns true if indices changed.
+        virtual bool valueWeights
+        (
+            const scalar t,
+            labelList& indices,
+            scalarField& weights
+        ) const = 0;
+
+        //- Calculate weights and indices to calculate integrand of t1..t2
+        //  from samples. Returns true if indices changed.
+        virtual bool integrationWeights
+        (
+            const scalar t1,
+            const scalar t2,
+            labelList& indices,
+            scalarField& weights
+        ) const = 0;
+
+        //- Helper: weighted sum
+        template<class ListType1, class ListType2>
+        static typename outerProduct
+        <
+            typename ListType1::value_type,
+            typename ListType2::value_type
+        >::type
+        weightedSum(const ListType1& f1, const ListType2& f2);
+
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "interpolationWeightsTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "interpolationWeights.H"
+#include "ListOps.H"
+#include "IOobject.H"
+#include "HashSet.H"
+#include "objectRegistry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ListType1, class ListType2>
+typename Foam::outerProduct
+<
+    typename ListType1::value_type,
+    typename ListType2::value_type
+>::type
+Foam::interpolationWeights::weightedSum
+(
+    const ListType1& f1,
+    const ListType2& f2
+)
+{
+    typedef typename outerProduct
+    <
+        typename ListType1::value_type,
+        typename ListType2::value_type
+    >::type returnType;
+
+    if (f1.size())
+    {
+        returnType SumProd = f1[0]*f2[0];
+        for (label i = 1; i < f1.size(); ++i)
+        {
+            SumProd += f1[i]*f2[i];
+        }
+        return SumProd;
+    }
+    else
+    {
+        return pTraits<returnType>::zero;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C

@@ -0,0 +1,249 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "linearInterpolationWeights.H"
+#include "addToRunTimeSelectionTable.H"
+#include "ListOps.H"
+#include "Pair.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(linearInterpolationWeights, 0);
+addToRunTimeSelectionTable
+(
+    interpolationWeights,
+    linearInterpolationWeights,
+    word
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::Pair<Foam::scalar> linearInterpolationWeights::integrationWeights
+(
+    const label i,
+    const scalar t
+) const
+{
+    // t is in range samples_[i] .. samples_[i+1]
+
+    scalar s = (t-samples_[i])/(samples_[i+1]-samples_[i]);
+
+    if (s < -SMALL || s > 1+SMALL)
+    {
+        FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
+            << "Value " << t << " outside range " << samples_[i]
+            << " .. " << samples_[i+1]
+            << exit(FatalError);
+    }
+
+    scalar d = samples_[i+1]-t;
+
+    return Pair<scalar>(d*0.5*(1-s), d*0.5*(1+s));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+linearInterpolationWeights::linearInterpolationWeights
+(
+    const scalarField& samples
+)
+:
+    interpolationWeights(samples),
+    index_(-1)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool linearInterpolationWeights::valueWeights
+(
+    const scalar t,
+    labelList& indices,
+    scalarField& weights
+) const
+{
+    bool indexChanged = false;
+
+    // Check if current timeIndex is still valid
+    if
+    (
+        index_ >= 0
+     && index_ < samples_.size()
+     && (
+            samples_[index_] <= t
+         && (index_ == samples_.size()-1 || t <= samples_[index_+1])
+        )
+    )
+    {
+        // index_ still at correct slot
+    }
+    else
+    {
+        // search for correct index
+        index_ = findLower(samples_, t);
+        indexChanged = true;
+    }
+
+
+    if (index_ == -1)
+    {
+        // Use first element only
+        indices.setSize(1);
+        weights.setSize(1);
+        indices[0] = 0;
+        weights[0] = 1.0;
+    }
+    else if (index_ == samples_.size()-1)
+    {
+        // Use last element only
+        indices.setSize(1);
+        weights.setSize(1);
+        indices[0] = samples_.size()-1;
+        weights[0] = 1.0;
+    }
+    else
+    {
+        // Interpolate
+        indices.setSize(2);
+        weights.setSize(2);
+
+        indices[0] = index_;
+        indices[1] = index_+1;
+
+        scalar t0 = samples_[indices[0]];
+        scalar t1 = samples_[indices[1]];
+        scalar deltaT = t1-t0;
+
+        weights[0] = (t1-t)/deltaT;
+        weights[1] = 1.0-weights[0];
+    }
+
+    return indexChanged;
+}
+
+
+bool linearInterpolationWeights::integrationWeights
+(
+    const scalar t1,
+    const scalar t2,
+    labelList& indices,
+    scalarField& weights
+) const
+{
+    if (t2 < t1-VSMALL)
+    {
+        FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
+            << "Integration should be in positive direction."
+            <<  " t1:" << t1 << " t2:" << t2
+            << exit(FatalError);
+    }
+
+    // Currently no fancy logic on cached index
+    label i1 = findLower(samples_, t1);
+    label i2 = findLower(samples_, t2);
+
+    // For now just fail if any outside table
+    if (i1 == -1 || i2 == samples_.size()-1)
+    {
+        FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
+            << "Integrating outside table " << samples_[0] << ".."
+            << samples_.last() << " not implemented."
+            << " t1:" << t1 << " t2:" << t2 << exit(FatalError);
+    }
+
+    label nIndices = i2-i1+2;
+
+
+    // Determine if indices already correct
+    bool anyChanged = false;
+
+    if (nIndices != indices.size())
+    {
+        anyChanged = true;
+    }
+    else
+    {
+        // Closer check
+
+        label index = i1;
+        forAll(indices, i)
+        {
+            if (indices[i] != index)
+            {
+                anyChanged = true;
+                break;
+            }
+            index++;
+        }
+    }
+
+    indices.setSize(nIndices);
+    weights.setSize(nIndices);
+    weights = 0.0;
+
+    // Sum from i1+1 to i2+1
+    for (label i = i1+1; i <= i2; i++)
+    {
+        scalar d = samples_[i+1]-samples_[i];
+        indices[i-i1] = i;
+        weights[i-i1] += 0.5*d;
+        indices[i+1-i1] = i+1;
+        weights[i+1-i1] += 0.5*d;
+    }
+
+    // Add from i1 to t1
+    {
+        Pair<scalar> i1Tot1 = integrationWeights(i1, t1);
+        indices[0] = i1;
+        weights[0] += i1Tot1.first();
+        indices[1] = i1+1;
+        weights[1] += i1Tot1.second();
+    }
+
+    // Subtract from t2 to i2+1
+    {
+        Pair<scalar> wghts = integrationWeights(i2, t2);
+        indices[i2-i1] = i2;
+        weights[i2-i1] += -wghts.first();
+        indices[i2-i1+1] = i2+1;
+        weights[i2-i1+1] += -wghts.second();
+    }
+
+    return anyChanged;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H

@@ -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/>.
+
+Class
+    Foam::linearInterpolationWeights
+
+Description
+
+SourceFiles
+    linearInterpolationWeights.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linearInterpolationWeights_H
+#define linearInterpolationWeights_H
+
+#include "interpolationWeights.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class linearInterpolationWeights Declaration
+\*---------------------------------------------------------------------------*/
+
+class linearInterpolationWeights
+:
+    public interpolationWeights
+{
+
+private:
+
+    // Private data
+
+        //- Cached index in samples from previous invocation
+        mutable label index_;
+
+    // Private Member Functions
+
+        //- Get weights of i and i+1 to calculate integration from t to
+        //  samples_[i+1]
+        Pair<scalar> integrationWeights
+        (
+            const label i,
+            const scalar t
+        ) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("linear");
+
+    // Constructors
+
+        //- Construct from components
+        linearInterpolationWeights
+        (
+            const scalarField& samples
+        );
+
+
+    //- Destructor
+    virtual ~linearInterpolationWeights()
+    {}
+
+
+    // Member Functions
+
+        //- Calculate weights and indices to calculate t from samples.
+        //  Returns true if indices changed.
+        virtual bool valueWeights
+        (
+            const scalar t,
+            labelList& indices,
+            scalarField& weights
+        ) const;
+
+        //- Calculate weights and indices to calculate integrand of t1..t2
+        //  from samples. Returns true if indices changed.
+        virtual bool integrationWeights
+        (
+            const scalar t1,
+            const scalar t2,
+            labelList& indices,
+            scalarField& weights
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C

@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "splineInterpolationWeights.H"
+#include "addToRunTimeSelectionTable.H"
+#include "ListOps.H"
+#include "linearInterpolationWeights.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(splineInterpolationWeights, 0);
+addToRunTimeSelectionTable
+(
+    interpolationWeights,
+    splineInterpolationWeights,
+    word
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+splineInterpolationWeights::splineInterpolationWeights
+(
+    const scalarField& samples,
+    const bool checkEqualDistance
+)
+:
+    interpolationWeights(samples),
+    index_(-1)
+{
+    if (checkEqualDistance && samples_.size() > 2)
+    {
+        const scalar interval = samples_[1]-samples[0];
+        for (label i = 2; i < samples_.size(); i++)
+        {
+            scalar d = samples_[i]-samples[i-1];
+
+            if (mag(d-interval) > SMALL)
+            {
+                WarningIn
+                (
+                    "splineInterpolationWeights::splineInterpolationWeights"
+                    "(const scalarField&)"
+                )   << "Spline interpolation only valid for constant intervals."
+                    << nl
+                    << "Interval 0-1 : " << interval << nl
+                    << "Interval " << i-1 << '-' << i << " : "
+                    << d << endl;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool splineInterpolationWeights::valueWeights
+(
+    const scalar t,
+    labelList& indices,
+    scalarField& weights
+) const
+{
+    bool indexChanged = false;
+
+    // linear interpolation
+    if (samples_.size() <= 2)
+    {
+        return linearInterpolationWeights(samples_).valueWeights
+        (
+            t,
+            indices,
+            weights
+        );
+    }
+
+    // Check if current timeIndex is still valid
+    if
+    (
+        index_ >= 0
+     && index_ < samples_.size()
+     && (
+            samples_[index_] <= t
+         && (index_ == samples_.size()-1 || t <= samples_[index_+1])
+        )
+    )
+    {
+        // index_ still at correct slot
+    }
+    else
+    {
+        // search for correct index
+        index_ = findLower(samples_, t);
+        indexChanged = true;
+    }
+
+
+    // Clamp if outside table
+    if (index_ == -1)
+    {
+        indices.setSize(1);
+        weights.setSize(1);
+
+        indices[0] = 0;
+        weights[0] = 1;
+        return indexChanged;
+    }
+    else if (index_ == samples_.size()-1)
+    {
+        indices.setSize(1);
+        weights.setSize(1);
+
+        indices[0] = samples_.size()-1;
+        weights[0] = 1;
+        return indexChanged;
+    }
+
+
+
+    label lo = index_;
+    label hi = index_+1;
+
+    // weighting
+    scalar mu = (t - samples_[lo])/(samples_[hi] - samples_[lo]);
+
+    scalar w0 = 0.5*(mu*(-1+mu*(2-mu)));            // coeff of lo-1
+    scalar w1 = 0.5*(2+mu*(mu*(-5 + mu*(3))));      // coeff of lo
+    scalar w2 = 0.5*(mu*(1 + mu*(4 + mu*(-3))));    // coeff of hi
+    scalar w3 = 0.5*(mu*mu*(-1 + mu));              // coeff of hi+1
+
+    if (lo > 0)
+    {
+        if (hi < samples_.size()-1)
+        {
+            // Four points available
+            indices.setSize(4);
+            weights.setSize(4);
+
+            indices[0] = lo-1;
+            indices[1] = lo;
+            indices[2] = hi;
+            indices[3] = hi+1;
+
+            weights[0] = w0;
+            weights[1] = w1;
+            weights[2] = w2;
+            weights[3] = w3;
+        }
+        else
+        {
+            // No y3 available. Extrapolate: y3=3*y2-y1
+            indices.setSize(3);
+            weights.setSize(3);
+
+            indices[0] = lo-1;
+            indices[1] = lo;
+            indices[2] = hi;
+
+            weights[0] = w0;
+            weights[1] = w1 - w3;
+            weights[2] = w2 + 2*w3;
+        }
+    }
+    else
+    {
+        // No y0 available. Extrapolate: y0=2*y1-y2;
+        if (hi < samples_.size()-1)
+        {
+            indices.setSize(3);
+            weights.setSize(3);
+
+            indices[0] = lo;
+            indices[1] = hi;
+            indices[2] = hi+1;
+
+            weights[0] = w1 + 2*w0;
+            weights[1] = w2 - w0;
+            weights[2] = w3;
+        }
+        else
+        {
+            indices.setSize(2);
+            weights.setSize(2);
+
+            indices[0] = lo;
+            indices[1] = hi;
+
+            weights[0] = w1 + 2*w0 - w3;
+            weights[1] = w2 - w0 + 2*w3;
+        }
+    }
+
+    return indexChanged;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::splineInterpolationWeights
+
+Description
+    Catmull-Rom spline interpolation.
+
+SourceFiles
+    splineInterpolationWeights.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef splineInterpolationWeights_H
+#define splineInterpolationWeights_H
+
+#include "interpolationWeights.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class splineInterpolationWeights Declaration
+\*---------------------------------------------------------------------------*/
+
+class splineInterpolationWeights
+:
+    public interpolationWeights
+{
+
+private:
+
+    // Private data
+
+        //- Cached index in samples from previous invocation
+        mutable label index_;
+
+public:
+
+    //- Runtime type information
+    TypeName("spline");
+
+    // Constructors
+
+        //- Construct from components. By default make sure samples are
+        //  equidistant.
+        splineInterpolationWeights
+        (
+            const scalarField& samples,
+            const bool checkEqualDistance = true
+        );
+
+
+    //- Destructor
+    virtual ~splineInterpolationWeights()
+    {}
+
+
+    // Member Functions
+
+        //- Calculate weights and indices to calculate t from samples.
+        //  Returns true if indices changed.
+        virtual bool valueWeights
+        (
+            const scalar t,
+            labelList& indices,
+            scalarField& weights
+        ) const;
+
+        //- Calculate weights and indices to calculate integrand of t1..t2
+        //  from samples. Returns true if indices changed.
+        virtual bool integrationWeights
+        (
+            const scalar t1,
+            const scalar t2,
+            labelList& indices,
+            scalarField& weights
+        ) const
+        {
+            notImplemented
+            (
+                "splineInterpolationWeights::integrationWeights(..)"
+            );
+            return false;
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/meshes/meshTools/mergePoints.C

@@ -180,7 +180,7 @@ bool Foam::mergePoints
     const bool verbose,
     labelList& pointMap,
     List<Type>& newPoints,
-    const Type& origin = Type::zero
+    const Type& origin
 )
 {
     label nUnique = mergePoints

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -44,7 +44,7 @@ License
 defineTypeNameAndDebug(Foam::globalMeshData, 0);
 
 // Geometric matching tolerance. Factor of mesh bounding box.
-const Foam::scalar Foam::globalMeshData::matchTol_ = 1E-8;
+const Foam::scalar Foam::globalMeshData::matchTol_ = 1e-8;
 
 namespace Foam
 {

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -893,7 +893,7 @@ Foam::polyMesh::cellTree() const
 
         Random rndGen(261782);
 
-        overallBb = overallBb.extend(rndGen, 1E-4);
+        overallBb = overallBb.extend(rndGen, 1e-4);
         overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
         overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 

src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.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
@@ -34,7 +34,7 @@ namespace Foam
 {
     defineTypeNameAndDebug(coupledPolyPatch, 0);
 
-    const scalar coupledPolyPatch::defaultMatchTol_ = 1E-4;
+    const scalar coupledPolyPatch::defaultMatchTol_ = 1e-4;
 
     template<>
     const char* NamedEnum<coupledPolyPatch::transformType, 4>::names[] =

src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C

@@ -26,6 +26,29 @@ License
 #include "TableBase.H"
 #include "Time.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+const Foam::interpolationWeights& Foam::TableBase<Type>::interpolator() const
+{
+    if (interpolatorPtr_.empty())
+    {
+        // Re-work table into linear list
+        tableSamples_.setSize(table_.size());
+        forAll(table_, i)
+        {
+            tableSamples_[i] = table_[i].first();
+        }
+        interpolatorPtr_ = interpolationWeights::New
+        (
+            interpolationScheme_,
+            tableSamples_
+        );
+    }
+    return interpolatorPtr_();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
@@ -39,6 +62,10 @@ Foam::TableBase<Type>::TableBase(const word& name, const dictionary& dict)
             dict.lookupOrDefault<word>("outOfBounds", "clamp")
         )
     ),
+    interpolationScheme_
+    (
+        dict.lookupOrDefault<word>("interpolationScheme", "linear")
+    ),
     table_(),
     dimensions_(dimless)
 {}
@@ -49,8 +76,11 @@ Foam::TableBase<Type>::TableBase(const TableBase<Type>& tbl)
 :
     name_(tbl.name_),
     boundsHandling_(tbl.boundsHandling_),
+    interpolationScheme_(tbl.interpolationScheme_),
     table_(tbl.table_),
-    dimensions_(tbl.dimensions_)
+    dimensions_(tbl.dimensions_),
+    tableSamples_(tbl.tableSamples_),
+    interpolatorPtr_(tbl.interpolatorPtr_)
 {}
 
 
@@ -307,6 +337,9 @@ void Foam::TableBase<Type>::convertTimeBase(const Time& t)
         scalar value = table_[i].first();
         table_[i].first() = t.userTimeToTime(value);
     }
+
+    tableSamples_.clear();
+    interpolatorPtr_.clear();
 }
 
 
@@ -325,88 +358,104 @@ Type Foam::TableBase<Type>::value(const scalar x) const
         return table_.last().second();
     }
 
-    // Find i such that x(i) < xDash < x(i+1)
+    // Use interpolator
-    label i = 0;
+    interpolator().valueWeights(x, currentIndices_, currentWeights_);
-    while ((table_[i+1].first() < xDash) && (i+1 < table_.size()))
+
+    Type t = currentWeights_[0]*table_[currentIndices_[0]].second();
+    for (label i = 1; i < currentIndices_.size(); i++)
     {
-        i++;
+        t += currentWeights_[i]*table_[currentIndices_[i]].second();
     }
+    return t;
 
-Info <<
+    //// Find i such that x(i) < xDash < x(i+1)
-     (xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
+    //label i = 0;
-      * (table_[i+1].second() - table_[i].second())
+    //while ((table_[i+1].first() < xDash) && (i+1 < table_.size()))
-      + table_[i].second() << endl;
+    //{
-
+    //    i++;
-    // Linear interpolation to find value
+    //}
-    return Type
+    //
-    (
+    //// Linear interpolation to find value
-        (xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
+    //return Type
-      * (table_[i+1].second() - table_[i].second())
+    //(
-      + table_[i].second()
+    //    (xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
-    );
+    //  * (table_[i+1].second() - table_[i].second())
+    //  + table_[i].second()
+    //);
 }
 
 
 template<class Type>
 Type Foam::TableBase<Type>::integrate(const scalar x1, const scalar x2) const
 {
-    // Initialise return value
+    // Use interpolator
-    Type sum = pTraits<Type>::zero;
+    interpolator().integrationWeights(x1, x2, currentIndices_, currentWeights_);
 
-    // Return zero if out of bounds
+    Type sum = currentWeights_[0]*table_[currentIndices_[0]].second();
-    if ((x1 > table_.last().first()) || (x2 < table_[0].first()))
+    for (label i = 1; i < currentIndices_.size(); i++)
     {
-        return sum;
+       sum += currentWeights_[i]*table_[currentIndices_[i]].second();
     }
-
-    // Find next index greater than x1
-    label id1 = 0;
-    while ((table_[id1].first() < x1) && (id1 < table_.size()))
-    {
-        id1++;
-    }
-
-    // Find next index less than x2
-    label id2 = table_.size() - 1;
-    while ((table_[id2].first() > x2) && (id2 >= 1))
-    {
-        id2--;
-    }
-
-    if ((id1 - id2) == 1)
-    {
-        // x1 and x2 lie within 1 interval
-        sum = 0.5*(value(x1) + value(x2))*(x2 - x1);
-    }
-    else
-    {
-        // x1 and x2 cross multiple intervals
-
-        // Integrate table body
-        for (label i=id1; i<id2; i++)
-        {
-            sum +=
-                (table_[i].second() + table_[i+1].second())
-              * (table_[i+1].first() - table_[i].first());
-        }
-        sum *= 0.5;
-
-        // Add table ends (partial segments)
-        if (id1 > 0)
-        {
-            sum += 0.5
-              * (value(x1) + table_[id1].second())
-              * (table_[id1].first() - x1);
-        }
-        if (id2 < table_.size() - 1)
-        {
-            sum += 0.5
-              * (table_[id2].second() + value(x2))
-              * (x2 - table_[id2].first());
-        }
-    }
-
     return sum;
+
+
+    //// Initialise return value
+    //Type sum = pTraits<Type>::zero;
+    //
+    //// Return zero if out of bounds
+    //if ((x1 > table_.last().first()) || (x2 < table_[0].first()))
+    //{
+    //    return sum;
+    //}
+    //
+    //// Find next index greater than x1
+    //label id1 = 0;
+    //while ((table_[id1].first() < x1) && (id1 < table_.size()))
+    //{
+    //    id1++;
+    //}
+    //
+    //// Find next index less than x2
+    //label id2 = table_.size() - 1;
+    //while ((table_[id2].first() > x2) && (id2 >= 1))
+    //{
+    //    id2--;
+    //}
+    //
+    //if ((id1 - id2) == 1)
+    //{
+    //    // x1 and x2 lie within 1 interval
+    //    sum = 0.5*(value(x1) + value(x2))*(x2 - x1);
+    //}
+    //else
+    //{
+    //    // x1 and x2 cross multiple intervals
+    //
+    //    // Integrate table body
+    //    for (label i=id1; i<id2; i++)
+    //    {
+    //        sum +=
+    //            (table_[i].second() + table_[i+1].second())
+    //          * (table_[i+1].first() - table_[i].first());
+    //    }
+    //    sum *= 0.5;
+    //
+    //    // Add table ends (partial segments)
+    //    if (id1 > 0)
+    //    {
+    //        sum += 0.5
+    //          * (value(x1) + table_[id1].second())
+    //          * (table_[id1].first() - x1);
+    //    }
+    //    if (id2 < table_.size() - 1)
+    //    {
+    //        sum += 0.5
+    //          * (table_[id2].second() + value(x2))
+    //          * (x2 - table_[id2].first());
+    //    }
+    //}
+    //
+    //return sum;
 }
 
 

src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H

@@ -38,6 +38,7 @@ SourceFiles
 #include "DataEntry.H"
 #include "Tuple2.H"
 #include "dimensionSet.H"
+#include "interpolationWeights.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,10 +81,13 @@ protected:
     // Protected data
 
         //- Table name
-        word name_;
+        const word name_;
 
         //- Enumeration for handling out-of-bound values
-        boundsHandling boundsHandling_;
+        const boundsHandling boundsHandling_;
+
+        //- Interpolation type
+        const word interpolationScheme_;
 
         //- Table data
         List<Tuple2<scalar, Type> > table_;
@@ -91,9 +95,24 @@ protected:
         //- The dimension set
         dimensionSet dimensions_;
 
+        //- Extracted values
+        mutable scalarField tableSamples_;
+
+        //- Interpolator method
+        mutable autoPtr<interpolationWeights> interpolatorPtr_;
+
+        //- Cached indices and weights
+        mutable labelList currentIndices_;
+
+        mutable scalarField currentWeights_;
+
 
     // Protected Member Functions
 
+
+        //- Return (demand driven) interpolator
+        const interpolationWeights& interpolator() const;
+
         //- Disallow default bitwise assignment
         void operator=(const TableBase<Type>&);
 

src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.C

@@ -45,6 +45,15 @@ Foam::TableFile<Type>::TableFile(const word& entryName, const dictionary& dict)
     fileName expandedFile(fName_);
     IFstream is(expandedFile.expand());
 
+    if (!is.good())
+    {
+        FatalIOErrorIn
+        (
+            "TableFile<Type>::TableFile(const word&, const dictionary&)",
+            is
+        )   << "Cannot open file." << exit(FatalIOError);
+    }
+
     is  >> this->table_;
 
     TableBase<Type>::check();

src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H

@@ -31,9 +31,10 @@ Description
         <entryName>   tableFile;
         tableFileCoeffs
         {
-            dimensions      [0 0 1 0 0]; // optional dimensions
+            dimensions          [0 0 1 0 0]; // optional dimensions
-            fileName        dataFile;    // name of data file
+            fileName            dataFile;    // name of data file
-            outOfBounds     clamp;       // optional out-of-bounds handling
+            outOfBounds         clamp;       // optional out-of-bounds handling
+            interpolationScheme linear;      // optional interpolation method
         }
     \endverbatim
 

src/combustionModels/FSD/FSD.H

@@ -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
@@ -48,7 +48,7 @@ Description
     release.
 
     If the turbulent fluctuation of the mixture fraction at the sub-grid level
-    is large (>1E-04) then a beta pdf is used for filtering.
+    is large (>1e-04) then a beta pdf is used for filtering.
 
     At the moment the flame area combustion model is only fit to work in a LES
     frame work. In RAS the subgrid fluctiuation has to be solved by an extra

src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C

@@ -240,10 +240,15 @@ Foam::dynamicRefineFvMesh::refine
         }
     }
 
+//    // Remove the stored tet base points
+//    tetBasePtIsPtr_.clear();
+//    // Remove the cell tree
+//    cellTreePtr_.clear();
 
     // Update fields
     updateMesh(map);
 
+
     // Move mesh
     /*
     pointField newPoints;

src/dynamicMesh/boundaryMesh/boundaryMesh.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,7 @@ defineTypeNameAndDebug(Foam::boundaryMesh, 0);
 const Foam::vector Foam::boundaryMesh::splitNormal_(3, 2, 1);
 
 // Distance to face tolerance for getNearest
-const Foam::scalar Foam::boundaryMesh::distanceTol_ = 1E-2;
+const Foam::scalar Foam::boundaryMesh::distanceTol_ = 1e-2;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -871,11 +871,11 @@ Foam::labelList Foam::boundaryMesh::getNearest
     {
         scalar sign = mesh().faceNormals()[bFaceI] & splitNormal_;
 
-        if (sign > -1E-5)
+        if (sign > -1e-5)
         {
             rightFaces.append(bFaceI);
         }
-        if (sign < 1E-5)
+        if (sign < 1e-5)
         {
             leftFaces.append(bFaceI);
         }
@@ -909,7 +909,7 @@ Foam::labelList Foam::boundaryMesh::getNearest
 
     // Extend domain slightly (also makes it 3D if was 2D)
     // Note asymmetry to avoid having faces align with octree cubes.
-    scalar tol = 1E-6 * overallBb.avgDim();
+    scalar tol = 1e-6 * overallBb.avgDim();
 
     point& bbMin = overallBb.min();
     bbMin.x() -= tol;

src/dynamicMesh/meshCut/cellLooper/geomCellLooper.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
@@ -39,7 +39,7 @@ License
 
 // Extension factor of edges to make sure we catch intersections through
 // edge endpoints
-const Foam::scalar Foam::geomCellLooper::pointEqualTol_ = 1E-3;
+const Foam::scalar Foam::geomCellLooper::pointEqualTol_ = 1e-3;
 
 
 // Snap cuts through edges onto edge endpoints. Fraction of edge length.

src/dynamicMesh/meshCut/directions/directions.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
@@ -125,7 +125,7 @@ void Foam::directions::check2D
 {
     if (correct2DPtr)
     {
-        if (mag(correct2DPtr->planeNormal() & vec) > 1E-6)
+        if (mag(correct2DPtr->planeNormal() & vec) > 1e-6)
         {
             FatalErrorIn("check2D") << "Specified vector " << vec
                 << "is not normal to plane defined in dynamicMeshDict."

src/dynamicMesh/motionSmoother/motionSmoother.C

@@ -900,7 +900,7 @@ Foam::tmp<Foam::scalarField> Foam::motionSmoother::movePoints
     {
         Pout<< "motionSmoother::movePoints : testing sync of newPoints."
             << endl;
-        testSyncPositions(newPoints, 1E-6*mesh_.bounds().mag());
+        testSyncPositions(newPoints, 1e-6*mesh_.bounds().mag());
     }
 
     // Move actual mesh points. Make sure to delete tetBasePtIs so it
@@ -1051,7 +1051,7 @@ bool Foam::motionSmoother::scaleMesh
                 totalDisplacement,
                 maxMagEqOp(),
                 vector::zero,   // null value
-                1E-6*mesh_.bounds().mag()
+                1e-6*mesh_.bounds().mag()
             );
         }
 

src/dynamicMesh/perfectInterface/perfectInterface.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
@@ -54,7 +54,7 @@ namespace Foam
 
 
 // Tolerance used as fraction of minimum edge length.
-const Foam::scalar Foam::perfectInterface::tol_ = 1E-3;
+const Foam::scalar Foam::perfectInterface::tol_ = 1e-3;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

src/dynamicMesh/polyMeshAdder/faceCoupleInfo.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
@@ -36,7 +36,7 @@ License
 
 defineTypeNameAndDebug(Foam::faceCoupleInfo, 0);
 
-const Foam::scalar Foam::faceCoupleInfo::angleTol_ = 1E-3;
+const Foam::scalar Foam::faceCoupleInfo::angleTol_ = 1e-3;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -1014,7 +1014,7 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
             mesh0,
             bndFaces                    // boundary faces only
         ),
-        overallBb.extend(rndGen, 1E-4), // overall search domain
+        overallBb.extend(rndGen, 1e-4), // overall search domain
         8,                              // maxLevel
         10,                             // leafsize
         3.0                             // duplicity

src/dynamicMesh/polyTopoChange/polyTopoChange/faceCollapser.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
@@ -277,7 +277,7 @@ void Foam::faceCollapser::setRefinement
             if (w <= dist[fpMin1])
             {
                 // Offset.
-                w = dist[fpMin1] + 1E-6*(dist[fpB] - dist[fpA]);
+                w = dist[fpMin1] + 1e-6*(dist[fpB] - dist[fpA]);
 
                 point newPoint
                 (
@@ -330,7 +330,7 @@ void Foam::faceCollapser::setRefinement
             if (w <= dist[fpMin1])
             {
                 // Offset.
-                w = dist[fpMin1] + 1E-6*(dist[fpB] - dist[fpA]);
+                w = dist[fpMin1] + 1e-6*(dist[fpB] - dist[fpA]);
 
                 point newPoint
                 (

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -4452,7 +4452,7 @@ void Foam::hexRef8::distribute(const mapDistributePolyMesh& map)
 
 void Foam::hexRef8::checkMesh() const
 {
-    const scalar smallDim = 1E-6 * mesh_.bounds().mag();
+    const scalar smallDim = 1e-6 * mesh_.bounds().mag();
 
     if (debug)
     {