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

applications/solvers/incompressible/pimpleFoam/pimpleFoam.C

@@ -30,7 +30,7 @@ Description
 
     Sub-models include:
     - turbulence modelling, i.e. laminar, RAS or LES
-    - run-time selectable fvOptions, e.g. MRF, explicit porosity
+    - run-time selectable finitie volume options, e.g. MRF, explicit porosity
 
 \*---------------------------------------------------------------------------*/
 

applications/solvers/lagrangian/reactingParcelFoam/Allwmake

@@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake
+wmake simpleReactingParcelFoam
 wmake LTSReactingParcelFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C

@@ -28,7 +28,7 @@ Description
     Local time stepping (LTS) solver for steady, compressible, laminar or
     turbulent reacting and non-reacting flow with multiphase Lagrangian
     parcels and porous media, including run-time selectable finitite volume
-    options, e.g. fvOptions, constraints
+    options, e.g. sources, constraints
 
     Note: ddtPhiCorr not used here when porous zones are active
     - not well defined for porous calculations

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -22,12 +22,12 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    porousExplicitSourceReactingParcelFoam
+    reactingParcelFoam
 
 Description
     Transient PIMPLE solver for compressible, laminar or turbulent flow with
     reacting multiphase Lagrangian parcels, including run-time selectable
-    finite volume options, e.g. fvOptions, constraints
+    finite volume options, e.g. sources, constraints
 
 \*---------------------------------------------------------------------------*/
 

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/EEqn.H

@@ -0,0 +1,32 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        mvConvection->fvmDiv(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        parcels.Sh(he)
+      + radiation->Sh(thermo)
+      + combustion->Sh()
+      + fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/files

@@ -0,0 +1,3 @@
+simpleReactingParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/simpleReactingParcelFoam

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options

@@ -0,0 +1,53 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lspecie \
+    -lfluidThermophysicalModels \
+    -lliquidProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lthermophysicalFunctions \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lradiationModels \
+    -lODE \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lcombustionModels \
+    -lfvOptions \
+    -lsampling

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H

@@ -0,0 +1,17 @@
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho.dimensionedInternalField()*g
+      + parcels.SU(U)
+      + fvOptions(rho, U)
+    );
+
+    UEqn().relax();
+
+    fvOptions.constrain(UEqn());
+
+    solve(UEqn() == -fvc::grad(p));
+
+    fvOptions.correct(U);

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H

@@ -0,0 +1,53 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    combustion->correct();
+    dQ = combustion->dQ();
+    label inertIndex = -1;
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YEqn
+            (
+                mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+             ==
+                parcels.SYi(i, Yi)
+              + combustion->R(Yi)
+              + fvOptions(rho, Yi)
+            );
+
+            YEqn.relax();
+
+            fvOptions.constrain(YEqn);
+
+            YEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H

@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicReactingMultiphaseCloud parcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    slgThermo
+);

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H

@@ -0,0 +1,98 @@
+    Info<< "Creating combustion model\n" << endl;
+
+    autoPtr<combustionModels::rhoCombustionModel> combustion
+    (
+        combustionModels::rhoCombustionModel::New(mesh)
+    );
+
+    rhoReactionThermo& thermo = combustion->thermo();
+    thermo.validate(args.executable(), "h", "e");
+
+    SLGThermo slgThermo(mesh, thermo);
+
+    basicMultiComponentMixture& composition = thermo.composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    const word inertSpecie(thermo.lookup("inertSpecie"));
+
+    if (!composition.contains(inertSpecie))
+    {
+        FatalErrorIn(args.executable())
+            << "Specified inert specie '" << inertSpecie << "' not found in "
+            << "species list. Available species:" << composition.species()
+            << exit(FatalError);
+    }
+
+    volScalarField& p = thermo.p();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    dimensionedScalar rhoMax(simple.dict().lookup("rhoMax"));
+    dimensionedScalar rhoMin(simple.dict().lookup("rhoMin"));
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    // Set the turbulence into the combustion model
+    combustion->setTurbulence(turbulence());
+
+    Info<< "Creating multi-variate interpolation scheme\n" << endl;
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll(Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(thermo.he());
+
+    volScalarField dQ
+    (
+        IOobject
+        (
+            "dQ",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    );

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H

@@ -0,0 +1,59 @@
+{
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    volScalarField rAU(1.0/UEqn().A());
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
+
+    UEqn.clear();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rho*rAU, p)
+         ==
+            parcels.Srho()
+          + fvOptions(psi, p, rho.name())
+        );
+
+        fvOptions.constrain(pEqn);
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+
+    p.relax();
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    #include "compressibleContinuityErrs.H"
+
+    U = HbyA - rAU*fvc::grad(p);
+    U.correctBoundaryConditions();
+    fvOptions.correct(U);
+
+    rho = thermo.rho();
+    rho = max(rho, rhoMin);
+    rho = min(rho, rhoMax);
+
+    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+}

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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
+    simpleReactingParcelFoam
+
+Description
+    Steady state SIMPLE solver for compressible, laminar or turbulent flow with
+    reacting multiphase Lagrangian parcels, including run-time selectable
+    finite volume options, e.g. sources, constraints
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulenceModel.H"
+#include "basicReactingMultiphaseCloud.H"
+#include "rhoCombustionModel.H"
+#include "radiationModel.H"
+#include "IOporosityModelList.H"
+#include "fvIOoptionList.H"
+#include "SLGThermo.H"
+#include "simpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+
+    simpleControl simple(mesh);
+
+    #include "createFields.H"
+    #include "createRadiationModel.H"
+    #include "createClouds.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        parcels.evolve();
+
+        // --- Pressure-velocity SIMPLE corrector loop
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "EEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H

@@ -0,0 +1,49 @@
+{
+    volScalarField k1
+    (
+        "k1",
+        alpha1*(thermo1.alpha()/rho1 + sqr(Ct)*nut2*thermo1.CpByCpv()/Prt)
+    );
+
+    volScalarField k2
+    (
+        "k2",
+        alpha2*(thermo2.alpha()/rho2 + nut2*thermo2.CpByCpv()/Prt)
+    );
+
+    volScalarField& he1 = thermo1.he();
+    volScalarField& he2 = thermo2.he();
+
+    Info<< max(he1) << min(he1) << endl;
+
+    fvScalarMatrix he1Eqn
+    (
+        fvm::ddt(alpha1, he1)
+      + fvm::div(alphaPhi1, he1)
+      - fvm::laplacian(k1, he1)
+     ==
+        heatTransferCoeff*(he1/thermo1.Cp())/rho1
+      - fvm::Sp(heatTransferCoeff/thermo1.Cp()/rho1, he1)
+      + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))
+    );
+
+    fvScalarMatrix he2Eqn
+    (
+        fvm::ddt(alpha2, he2)
+      + fvm::div(alphaPhi2, he2)
+      - fvm::laplacian(k2, he2)
+     ==
+        heatTransferCoeff*(he2/thermo2.Cp())/rho2
+      - fvm::Sp(heatTransferCoeff/thermo2.Cp()/rho2, he2)
+      + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))
+    );
+
+    he1Eqn.relax();
+    he1Eqn.solve();
+
+    he2Eqn.relax();
+    he2Eqn.solve();
+
+    thermo1.correct();
+    thermo2.correct();
+}

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -IturbulenceModel \
@@ -8,6 +9,8 @@ EXE_INC = \
     -Iaveraging
 
 EXE_LIBS = \
+    -lfluidThermophysicalModels \
+    -lspecie \
     -lincompressibleTransportModels \
     -lcompressiblePhaseModel \
     -lcompressibleEulerianInterfacialModels \

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H

@@ -1,36 +0,0 @@
-{
-    volScalarField kByCp1("kByCp1", alpha1*(k1/Cp1/rho1 + sqr(Ct)*nut2/Prt));
-    volScalarField kByCp2("kByCp2", alpha2*(k2/Cp2/rho2 + nut2/Prt));
-
-    fvScalarMatrix T1Eqn
-    (
-        fvm::ddt(alpha1, T1)
-      + fvm::div(alphaPhi1, T1)
-      - fvm::laplacian(kByCp1, T1)
-     ==
-        heatTransferCoeff*T2/Cp1/rho1
-      - fvm::Sp(heatTransferCoeff/Cp1/rho1, T1)
-      + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))/Cp1
-    );
-
-    fvScalarMatrix T2Eqn
-    (
-        fvm::ddt(alpha2, T2)
-      + fvm::div(alphaPhi2, T2)
-      - fvm::laplacian(kByCp2, T2)
-     ==
-        heatTransferCoeff*T1/Cp2/rho2
-      - fvm::Sp(heatTransferCoeff/Cp2/rho2, T2)
-      + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))/Cp2
-    );
-
-    T1Eqn.relax();
-    T1Eqn.solve();
-
-    T2Eqn.relax();
-    T2Eqn.solve();
-
-    // Update compressibilities
-    psi1 = 1.0/(R1*T1);
-    psi2 = 1.0/(R2*T2);
-}

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H

@@ -12,12 +12,12 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         }
         else // If not using kinetic theory is using Ct model
         {
-            nuEff1 = sqr(Ct)*nut2 + nu1;
+            nuEff1 = sqr(Ct)*nut2 + thermo1.mu()/rho1;
         }
 
         volTensorField Rc1
         (
-            "Rc1",
+            "Rc",
             (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T
         );
 
@@ -52,7 +52,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         volTensorField gradU2T(fvc::grad(U2)().T());
         volTensorField Rc2
         (
-            "Rc2",
+            "Rc",
             (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T
         );
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H

@@ -1,9 +1,9 @@
-surfaceScalarField alphaPhi1("alphaPhi1", phi1);
-surfaceScalarField alphaPhi2("alphaPhi2", phi2);
+surfaceScalarField alphaPhi1("alphaPhi", phi1);
+surfaceScalarField alphaPhi2("alphaPhi", phi2);
 
 {
-    word scheme("div(phi,alpha1)");
-    word schemer("div(phir,alpha1)");
+    word scheme("div(phi,alpha)");
+    word schemer("div(phir,alpha)");
 
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
@@ -13,7 +13,7 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
         surfaceScalarField alpha1f(fvc::interpolate(alpha1));
         surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf());
         phir += phipp;
-        phic += fvc::interpolate(alpha1)*phipp;
+        phic += alpha1f*phipp;
     }
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -68,16 +68,22 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
 
         if (g0.value() > 0.0)
         {
+            surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+
             ppMagf =
                 fvc::interpolate((1.0/rho1)*rAU1)
-               *fvc::interpolate
-                (
-                    g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
-                );
+               *g0*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
+
+            // ppMagf =
+            //     fvc::interpolate((1.0/rho1)*rAU1)
+            //    *fvc::interpolate
+            //     (
+            //         g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
+            //     );
 
             alpha1Eqn -= fvm::laplacian
             (
-                (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
+                alpha1f*ppMagf,
                 alpha1,
                 "laplacian(alphaPpMag,alpha1)"
             );

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "rhoThermo.H"
 #include "nearWallDist.H"
 #include "wallFvPatch.H"
 #include "fixedValueFvsPatchFields.H"
@@ -86,7 +87,7 @@ int main(int argc, char *argv[])
             #include "alphaEqn.H"
             #include "kEpsilon.H"
             #include "interfacialCoeffs.H"
-            #include "TEqns.H"
+            #include "EEqns.H"
             #include "UEqns.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H

@@ -12,55 +12,43 @@
         )
     );
 
+    word phase1Name
+    (
+        transportProperties.found("phases")
+      ? wordList(transportProperties.lookup("phases"))[0]
+      : "phase1"
+    );
+
+    word phase2Name
+    (
+        transportProperties.found("phases")
+      ? wordList(transportProperties.lookup("phases"))[1]
+      : "phase2"
+    );
+
     autoPtr<phaseModel> phase1 = phaseModel::New
     (
         mesh,
         transportProperties,
-        "1"
+        phase1Name
     );
 
     autoPtr<phaseModel> phase2 = phaseModel::New
     (
         mesh,
         transportProperties,
-        "2"
+        phase2Name
     );
 
+    volScalarField& alpha1 = phase1();
+    volScalarField& alpha2 = phase2();
+    alpha2 = scalar(1) - alpha1;
+
     volVectorField& U1 = phase1->U();
     surfaceScalarField& phi1 = phase1->phi();
-    const dimensionedScalar& nu1 = phase1->nu();
-    const dimensionedScalar& k1 = phase1->kappa();
-    const dimensionedScalar& Cp1 = phase1->Cp();
-    dimensionedScalar rho10
-    (
-        "rho0",
-        dimDensity,
-        transportProperties.subDict("phase1").lookup("rho0")
-    );
-    dimensionedScalar R1
-    (
-        "R",
-        dimensionSet(0, 2, -2, -1, 0),
-        transportProperties.subDict("phase1").lookup("R")
-    );
 
     volVectorField& U2 = phase2->U();
     surfaceScalarField& phi2 = phase2->phi();
-    const dimensionedScalar& nu2 = phase2->nu();
-    const dimensionedScalar& k2 = phase2->kappa();
-    const dimensionedScalar& Cp2 = phase2->Cp();
-    dimensionedScalar rho20
-    (
-        "rho0",
-        dimDensity,
-        transportProperties.subDict("phase2").lookup("rho0")
-    );
-    dimensionedScalar R2
-    (
-        "R",
-        dimensionSet(0, 2, -2, -1, 0),
-        transportProperties.subDict("phase2").lookup("R")
-    );
 
     dimensionedScalar pMin
     (
@@ -69,102 +57,16 @@
         transportProperties.lookup("pMin")
     );
 
+    rhoThermo& thermo1 = phase1->thermo();
+    rhoThermo& thermo2 = phase2->thermo();
 
-    Info<< "Reading field T1" << endl;
-    volScalarField T1
-    (
-        IOobject
-        (
-            "T1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
+    volScalarField& p = thermo1.p();
 
-    Info<< "Reading field T2" << endl;
-    volScalarField T2
-    (
-        IOobject
-        (
-            "T2",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
+    volScalarField rho1("rho" + phase1Name, thermo1.rho());
+    const volScalarField& psi1 = thermo1.psi();
 
-    volScalarField psi1
-    (
-        IOobject
-        (
-            "psi1",
-            runTime.timeName(),
-            mesh
-        ),
-        1.0/(R1*T1)
-    );
-
-    volScalarField psi2
-    (
-        IOobject
-        (
-            "psi2",
-            runTime.timeName(),
-            mesh
-        ),
-        1.0/(R2*T2)
-    );
-    psi2.oldTime();
-
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField rho1("rho1", rho10 + psi1*p);
-    volScalarField rho2("rho2", rho20 + psi2*p);
-
-
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField alpha2
-    (
-        IOobject
-        (
-            "alpha2",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        scalar(1) - alpha1
-    );
+    volScalarField rho2("rho" + phase2Name, thermo2.rho());
+    const volScalarField& psi2 = thermo2.psi();
 
     volVectorField U
     (
@@ -179,27 +81,6 @@
         alpha1*U1 + alpha2*U2
     );
 
-    dimensionedScalar Cvm
-    (
-        "Cvm",
-        dimless,
-        transportProperties.lookup("Cvm")
-    );
-
-    dimensionedScalar Cl
-    (
-        "Cl",
-        dimless,
-        transportProperties.lookup("Cl")
-    );
-
-    dimensionedScalar Ct
-    (
-        "Ct",
-        dimless,
-        transportProperties.lookup("Ct")
-    );
-
     surfaceScalarField phi
     (
         IOobject
@@ -226,8 +107,6 @@
         alpha1*rho1 + alpha2*rho2
     );
 
-    #include "createRASTurbulence.H"
-
     Info<< "Calculating field DDtU1 and DDtU2\n" << endl;
 
     volVectorField DDtU1
@@ -248,6 +127,29 @@
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
 
+    dimensionedScalar Cvm
+    (
+        "Cvm",
+        dimless,
+        transportProperties.lookup("Cvm")
+    );
+
+    dimensionedScalar Cl
+    (
+        "Cl",
+        dimless,
+        transportProperties.lookup("Cl")
+    );
+
+    dimensionedScalar Ct
+    (
+        "Ct",
+        dimless,
+        transportProperties.lookup("Ct")
+    );
+
+    #include "createRASTurbulence.H"
+
     IOdictionary interfacialProperties
     (
         IOobject
@@ -297,8 +199,8 @@
     if
     (
         !(
-            dispersedPhase == "1"
-         || dispersedPhase == "2"
+            dispersedPhase == phase1Name
+         || dispersedPhase == phase2Name
          || dispersedPhase == "both"
         )
     )
@@ -337,7 +239,7 @@
     (
         IOobject
         (
-            "rAU1",
+            "rAU" + phase1Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -387,5 +289,5 @@
     );
 
     Info<< "Creating field kinetic energy K\n" << endl;
-    volScalarField K1("K1", 0.5*magSqr(U1));
-    volScalarField K2("K2", 0.5*magSqr(U2));
+    volScalarField K1("K" + phase1Name, 0.5*magSqr(U1));
+    volScalarField K2("K" + phase2Name, 0.5*magSqr(U2));

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H

@@ -151,7 +151,7 @@
     (
         IOobject
         (
-            "nut2",
+            "nut" + phase2Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -165,13 +165,13 @@
     (
         IOobject
         (
-            "nuEff1",
+            "nuEff" + phase1Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        sqr(Ct)*nut2 + nu1
+        sqr(Ct)*nut2 + thermo1.mu()/rho1
     );
 
     Info<< "Calculating field nuEff2\n" << endl;
@@ -179,11 +179,11 @@
     (
         IOobject
         (
-            "nuEff2",
+            "nuEff" + phase2Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        nut2 + nu2
+        nut2 + thermo2.mu()/rho2
     );

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H

@@ -38,12 +38,12 @@ volScalarField heatTransferCoeff
     volVectorField Ur(U1 - U2);
     volScalarField magUr(mag(Ur) + residualSlip);
 
-    if (dispersedPhase == "1")
+    if (dispersedPhase == phase1Name)
     {
         dragCoeff = drag1->K(magUr);
         heatTransferCoeff = heatTransfer1->K(magUr);
     }
-    else if (dispersedPhase == "2")
+    else if (dispersedPhase == phase2Name)
     {
         dragCoeff = drag2->K(magUr);
         heatTransferCoeff = heatTransfer2->K(magUr);

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options

@@ -1,7 +1,9 @@
 EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I../phaseModel/lnInclude
 
 LIB_LIBS = \
-    -lcompressiblePhaseModel
-
+    -lcompressiblePhaseModel \
+    -lfluidThermophysicalModels \
+    -lspecie

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,8 +72,7 @@ Foam::tmp<Foam::volScalarField> Foam::heatTransferModels::RanzMarshall::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    dimensionedScalar Prb =
-        phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa();
+    volScalarField Prb(phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa());
     volScalarField Nu(scalar(2) + 0.6*sqrt(Re)*cbrt(Prb));
 
     return 6.0*phase2_.kappa()*Nu/sqr(phase1_.d());

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/foam/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I../phaseModel/lnInclude \
     -I../interfacialModels/lnInclude

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,8 +45,8 @@ Foam::kineticTheoryModel::kineticTheoryModel
     phi1_(phase1.phi()),
     draga_(draga),
 
-    rho1_(phase1.rho()),
-    nu1_(phase1.nu()),
+    rho1_(phase1.rho()[0]), //***HGW
+    nu1_(phase1.nu()()[0]), //***HGW
 
     kineticTheoryProperties_
     (
@@ -120,7 +120,7 @@ Foam::kineticTheoryModel::kineticTheoryModel
     (
         IOobject
         (
-            "mu1",
+            "mu" + phase1.name(),
             U1_.time().timeName(),
             U1_.mesh(),
             IOobject::NO_READ,

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -1,6 +1,6 @@
 {
-    rho1 = rho10 + psi1*p;
-    rho2 = rho20 + psi2*p;
+    rho1 = thermo1.rho();
+    rho2 = thermo2.rho();
 
     surfaceScalarField alpha1f(fvc::interpolate(alpha1));
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
@@ -11,10 +11,10 @@
     surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rAU1));
     surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rAU2));
 
-    volVectorField HbyA1("HbyA1", U1);
+    volVectorField HbyA1("HbyA" + phase1Name, U1);
     HbyA1 = rAU1*U1Eqn.H();
 
-    volVectorField HbyA2("HbyA2", U2);
+    volVectorField HbyA2("HbyA" + phase2Name, U2);
     HbyA2 = rAU2*U2Eqn.H();
 
     mrfZones.absoluteFlux(phi1.oldTime());
@@ -38,14 +38,14 @@
 
     surfaceScalarField phiHbyA1
     (
-        "phiHbyA1",
+        "phiHbyA" + phase1Name,
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
     );
 
     surfaceScalarField phiHbyA2
     (
-        "phiHbyA2",
+        "phiHbyA" + phase2Name,
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
     );
@@ -96,8 +96,16 @@
     //}
     //else
     {
-        surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi1);
-        surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2);
+        surfaceScalarField phid1
+        (
+            "phid" + phase1Name,
+            fvc::interpolate(psi1)*phi1
+        );
+        surfaceScalarField phid2
+        (
+            "phid" + phase2Name,
+            fvc::interpolate(psi2)*phi2
+        );
 
         pEqnComp1 =
             fvc::ddt(rho1) + psi1*correction(fvm::ddt(p))
@@ -174,13 +182,15 @@
 
     p = max(p, pMin);
 
-    rho1 = rho10 + psi1*p;
-    rho2 = rho20 + psi2*p;
+    thermo1.correct();
+    thermo2.correct();
+    rho1 = thermo1.rho();
+    rho2 = thermo2.rho();
 
     K1 = 0.5*magSqr(U1);
     K2 = 0.5*magSqr(U2);
 
-    //***HGW if (thermo.dpdt())
+    if (thermo1.dpdt())
     {
         dpdt = fvc::ddt(p);
     }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options

@@ -1,6 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude
 
 LIB_LIBS = \
-    -lincompressibleTransportModels
+    -lincompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,35 +37,22 @@ Foam::phaseModel::phaseModel
     const word& phaseName
 )
 :
-    dict_
+    volScalarField
     (
-        transportProperties.subDict("phase" + phaseName)
+        IOobject
+        (
+            "alpha" + phaseName,
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("alpha", dimless, 0)
     ),
     name_(phaseName),
-    nu_
-    (
-        "nu",
-        dimensionSet(0, 2, -1, 0, 0),
-        dict_.lookup("nu")
-    ),
-    kappa_
-    (
-        "kappa",
-        dimensionSet(1, 1, -3, -1, 0),
-        dict_.lookup("kappa")
-    ),
-    Cp_
-    (
-        "Cp",
-        dimensionSet(0, 2, -2, -1, 0),
-        dict_.lookup("Cp")
-    ),
-    rho_
-    (
-        "rho",
-        dimDensity,
-        dict_.lookup("rho")
-    ),
+    phaseDict_(transportProperties.subDict(phaseName)),
+    thermo_(rhoThermo::New(mesh, phaseName)),
     U_
     (
         IOobject
@@ -79,6 +66,8 @@ Foam::phaseModel::phaseModel
         mesh
     )
 {
+    thermo_->validate("phaseModel " + phaseName, "h", "e");
+
     const word phiName = "phi" + phaseName;
 
     IOobject phiHeader
@@ -147,7 +136,7 @@ Foam::phaseModel::phaseModel
 
     dPtr_ = diameterModel::New
     (
-        dict_,
+        phaseDict_,
         *this
     );
 }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ SourceFiles
 #include "dimensionedScalar.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "rhoThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,25 +51,18 @@ class diameterModel;
 \*---------------------------------------------------------------------------*/
 
 class phaseModel
+:
+    public volScalarField
 {
     // Private data
 
-        dictionary dict_;
-
         //- Name of phase
         word name_;
 
-        //- Kinematic viscosity
-        dimensionedScalar nu_;
+        dictionary phaseDict_;
 
-        //- Thermal conductivity
-        dimensionedScalar kappa_;
-
-        //- Heat capacity
-        dimensionedScalar Cp_;
-
-        //- Density
-        dimensionedScalar rho_;
+        //- Thermophysical properties
+        autoPtr<rhoThermo> thermo_;
 
         //- Velocity
         volVectorField U_;
@@ -116,24 +110,34 @@ public:
 
         tmp<volScalarField> d() const;
 
-        const dimensionedScalar& nu() const
+        tmp<volScalarField> nu() const
         {
-            return nu_;
+            return thermo_->mu()/thermo_->rho();
         }
 
-        const dimensionedScalar& kappa() const
+        tmp<volScalarField> kappa() const
         {
-            return kappa_;
+            return thermo_->kappa();
         }
 
-        const dimensionedScalar& Cp() const
+        tmp<volScalarField> Cp() const
         {
-            return Cp_;
+            return thermo_->Cp();
         }
 
-        const dimensionedScalar& rho() const
+        const volScalarField& rho() const
         {
-            return rho_;
+            return thermo_->rho();
+        }
+
+        const rhoThermo& thermo() const
+        {
+            return thermo_();
+        }
+
+        rhoThermo& thermo()
+        {
+            return thermo_();
         }
 
         const volVectorField& U() const

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H

@@ -61,4 +61,4 @@ if (turbulence)
     #include "wallViscosity.H"
 }
 
-nuEff2 = nut2 + nu2;
+nuEff2 = nut2 + thermo2.mu()/rho2;

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H

@@ -4,7 +4,6 @@
     scalar Cmu25 = ::pow(Cmu.value(), 0.25);
     scalar Cmu75 = ::pow(Cmu.value(), 0.75);
     scalar kappa_ = kappa.value();
-    scalar nu2_ = nu2.value();
 
     const fvPatchList& patches = mesh.boundary();
 
@@ -30,6 +29,8 @@
     forAll(patches, patchi)
     {
         const fvPatch& currPatch = patches[patchi];
+        const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi];
+        const scalarField& rho2_ = rho2.boundaryField()[patchi];
 
         if (isA<wallFvPatch>(currPatch))
         {
@@ -52,7 +53,8 @@
                     /(kappa_*y[patchi][facei]);
 
                 G[faceCelli] +=
-                    (nut2w[facei] + nu2_)*magFaceGradU[facei]
+                    (nut2w[facei] + mu2_[facei]/rho2_[facei])
+                   *magFaceGradU[facei]
                    *Cmu25*::sqrt(k[faceCelli])
                    /(kappa_*y[patchi][facei]);
             }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H

@@ -2,13 +2,14 @@
     scalar Cmu25 = ::pow(Cmu.value(), 0.25);
     scalar kappa_ = kappa.value();
     scalar E_ = E.value();
-    scalar nu2_ = nu2.value();
 
     const fvPatchList& patches = mesh.boundary();
 
     forAll(patches, patchi)
     {
         const fvPatch& currPatch = patches[patchi];
+        const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi];
+        const scalarField& rho2_ = rho2.boundaryField()[patchi];
 
         if (isA<wallFvPatch>(currPatch))
         {
@@ -20,11 +21,14 @@
 
                 // calculate yPlus
                 scalar yPlus =
-                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nu2_;
+                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
+                   /(mu2_[facei]/rho2_[facei]);
 
                 if (yPlus > 11.6)
                 {
-                    nutw[facei] = nu2_*(yPlus*kappa_/::log(E_*yPlus) - 1);
+                    nutw[facei] =
+                        (mu2_[facei]/rho2_[facei])
+                       *(yPlus*kappa_/::log(E_*yPlus) - 1);
                 }
                 else
                 {

applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H

@@ -13,7 +13,7 @@
         surfaceScalarField alpha1f(fvc::interpolate(alpha1));
         surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf());
         phir += phipp;
-        phic += fvc::interpolate(alpha1)*phipp;
+        phic += alpha1f*phipp;
     }
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -52,18 +52,23 @@
 
         if (g0.value() > 0)
         {
-            ppMagf = rAU1f*fvc::interpolate
-            (
-                (1.0/(rho1*(alpha1 + scalar(0.0001))))
-               *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
-            );
+            surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+
+            // ppMagf = rAU1f*fvc::interpolate
+            // (
+            //     (1.0/(rho1*(alpha1 + scalar(0.0001))))
+            //    *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
+            // );
+            ppMagf =
+                rAU1f/(alpha1f + scalar(0.0001))
+               *(g0/rho1)*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
 
             fvScalarMatrix alpha1Eqn
             (
                 fvm::ddt(alpha1) - fvc::ddt(alpha1)
               - fvm::laplacian
                 (
-                    (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
+                    alpha1f*ppMagf,
                     alpha1,
                     "laplacian(alpha1PpMag,alpha1)"
                 )

applications/test/Circulator/Test-Circulator.C

@@ -38,94 +38,6 @@ Description
 using namespace Foam;
 
 
-// return
-//   0: no match
-//  +1: identical
-//  -1: same face, but different orientation
-label compare(const face& a, const face& b)
-{
-    // Basic rule: we assume that the sequence of labels in each list
-    // will be circular in the same order (but not necessarily in the
-    // same direction or from the same starting point).
-
-    // Trivial reject: faces are different size
-    label sizeA = a.size();
-    label sizeB = b.size();
-
-    if (sizeA != sizeB || sizeA == 0)
-    {
-        return 0;
-    }
-
-    const_circulator<face> aCirc(a);
-    const_circulator<face> bCirc(b);
-
-    // Rotate face b until its element matches the starting element of face a.
-    do
-    {
-        if (aCirc() == bCirc())
-        {
-            // Set bCirc fulcrum to its iterator and increment the iterators
-            bCirc.setFulcrumToIterator();
-            ++aCirc;
-            ++bCirc;
-
-            break;
-        }
-    } while (bCirc.circulate(CirculatorBase::CLOCKWISE));
-
-    // Look forwards around the faces for a match
-    do
-    {
-        if (aCirc() != bCirc())
-        {
-            break;
-        }
-    }
-    while
-    (
-        aCirc.circulate(CirculatorBase::CLOCKWISE),
-        bCirc.circulate(CirculatorBase::CLOCKWISE)
-    );
-
-    // If the circulator has stopped then faces a and b matched.
-    if (!aCirc.circulate())
-    {
-        return 1;
-    }
-    else
-    {
-        // Reset the circulators back to their fulcrum
-        aCirc.setIteratorToFulcrum();
-        bCirc.setIteratorToFulcrum();
-        ++aCirc;
-        --bCirc;
-    }
-
-    // Look backwards around the faces for a match
-    do
-    {
-        if (aCirc() != bCirc())
-        {
-            break;
-        }
-    }
-    while
-    (
-        aCirc.circulate(CirculatorBase::CLOCKWISE),
-        bCirc.circulate(CirculatorBase::ANTICLOCKWISE)
-    );
-
-    // If the circulator has stopped then faces a and b matched.
-    if (!aCirc.circulate())
-    {
-        return -1;
-    }
-
-    return 0;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
 
@@ -184,40 +96,44 @@ int main(int argc, char *argv[])
 
     Info<< nl << nl << "Compare two faces: " << endl;
     face a(identity(5));
-    Info<< "Compare " << a << " and " << a << " Match = " << compare(a, a)
+    Info<< "Compare " << a << " and " << a << " Match = " << face::compare(a, a)
         << endl;
 
     face b(reverseList(a));
-    Info<< "Compare " << a << " and " << b << " Match = " << compare(a, b)
+    Info<< "Compare " << a << " and " << b << " Match = " << face::compare(a, b)
         << endl;
 
     face c(a);
     c[4] = 3;
-    Info<< "Compare " << a << " and " << c << " Match = " << compare(a, c)
+    Info<< "Compare " << a << " and " << c << " Match = " << face::compare(a, c)
         << endl;
 
     face d(rotateList(a, 2));
-    Info<< "Compare " << a << " and " << d << " Match = " << compare(a, d)
+    Info<< "Compare " << a << " and " << d << " Match = " << face::compare(a, d)
         << endl;
 
     face g(labelList(5, 1));
     face h(g);
-    Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+    Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h)
         << endl;
 
     g[0] = 2;
     h[3] = 2;
-    Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+    Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h)
         << endl;
 
     g[4] = 3;
     h[4] = 3;
-    Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+    Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h)
         << endl;
 
     face face1(identity(1));
     Info<< "Compare " << face1 << " and " << face1
-        << " Match = " << compare(face1, face1) << endl;
+        << " Match = " << face::compare(face1, face1) << endl;
+
+    face face2(identity(1)+1);
+    Info<< "Compare " << face1 << " and " << face2
+        << " Match = " << face::compare(face1, face2) << endl;
 
     Info<< nl << nl << "Zero face" << nl << endl;
 

applications/utilities/mesh/advanced/collapseEdges/collapseDict

@@ -17,7 +17,7 @@ FoamFile
 collapseEdgesCoeffs
 {
     // Edges shorter than this absolute value will be merged
-    minimumEdgeLength   1e-8;
+    minimumEdgeLength   1e-6;
 
     // The maximum angle between two edges that share a point attached to
     // no other edges
@@ -25,7 +25,7 @@ collapseEdgesCoeffs
 
     // The amount that minimumEdgeLength will be reduced by for each
     // edge if that edge's collapse generates a poor quality face
-    reductionFactor                         0.5;
+    reductionFactor     0.5;
 }
 
 
@@ -67,14 +67,17 @@ meshQualityCoeffs
 {
     // Name of the dictionary that has the mesh quality coefficients used
     // by motionSmoother::checkMesh
-    meshQualityCoeffDict    meshQualityDict;
+    #include                    "meshQualityDict";
+
+    // Maximum number of smoothing iterations for the reductionFactors
+    maximumSmoothingIterations  2;
 
     // Maximum number of outer iterations is mesh quality checking is enabled
-    maximumIterations       10;
+    maximumIterations           10;
 
     // Maximum number of iterations deletion of a point can cause a bad face
     // to be constructed before it is forced to not be deleted
-    maxPointErrorCount      5;
+    maxPointErrorCount          5;
 }
 
 

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict

@@ -70,6 +70,9 @@ extrudeModel    linearNormal;
 //- Extrudes into sphere around (0 0 0)
 // extrudeModel   linearRadial;
 
+//- Extrudes into sphere around (0 0 0) with specified radii
+//extrudeModel        radial;
+
 //- Extrudes into sphere with grading according to pressure (atmospherics)
 // extrudeModel   sigmaRadial;
 
@@ -98,6 +101,14 @@ linearDirectionCoeffs
 linearRadialCoeffs
 {
      R              0.1;
+    // Optional inner radius
+    Rsurface        0.01;
+}
+
+radialCoeffs
+{
+    // Radii specified through interpolation table
+    R               table ((0 0.01)(3 0.03)(10 0.1));
 }
 
 sigmaRadialCoeffs
@@ -107,4 +118,5 @@ sigmaRadialCoeffs
     pStrat          1;
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -1,5 +1,4 @@
 EXE_INC = \
-    /* -DFULLDEBUG -g -O0 */ \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \

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

@@ -26,12 +26,6 @@ License
 #include "patchToPoly2DMesh.H"
 #include "PatchTools.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::patchToPoly2DMesh::flipFaceOrder()
@@ -275,12 +269,9 @@ void Foam::patchToPoly2DMesh::createPolyMeshComponents()
     addPatchFacesToOwner();
 }
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-//- Construct from a primitivePatch
 Foam::patchToPoly2DMesh::patchToPoly2DMesh
 (
     const MeshedSurface<face>& patch,
@@ -321,13 +312,5 @@ void Foam::patchToPoly2DMesh::createMesh()
     }
 }
 
-// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
-
 
 // ************************************************************************* //

applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H

@@ -51,6 +51,7 @@ class patchToPoly2DMesh
 {
     // Private data
 
+        // Reference to the meshed surface
         const MeshedSurface<face>& patch_;
 
         const wordList& patchNames_;
@@ -62,11 +63,16 @@ class patchToPoly2DMesh
         const EdgeMap<label>& mapEdgesRegion_;
 
         pointField points_;
+
         faceList faces_;
+
         labelList owner_;
+
         labelList neighbour_;
 
-        //- Description of data_
+
+    // Private Member Functions
+
         void flipFaceOrder();
 
         void createNeighbours();
@@ -79,9 +85,6 @@ class patchToPoly2DMesh
 
         void createPolyMeshComponents();
 
-
-    // Private Member Functions
-
         //- Disallow default bitwise copy construct
         patchToPoly2DMesh(const patchToPoly2DMesh&);
 
@@ -110,47 +113,47 @@ public:
     // Member Functions
 
         // Access
-        pointField& points()
-        {
-            return points_;
-        }
 
-        faceList& faces()
-        {
-            return faces_;
-        }
+            pointField& points()
+            {
+                return points_;
+            }
 
-        labelList& owner()
-        {
-            return owner_;
-        }
+            faceList& faces()
+            {
+                return faces_;
+            }
 
-        labelList& neighbour()
-        {
-            return neighbour_;
-        }
+            labelList& owner()
+            {
+                return owner_;
+            }
 
-        const wordList& patchNames() const
-        {
-            return patchNames_;
-        }
+            labelList& neighbour()
+            {
+                return neighbour_;
+            }
 
-        const labelList& patchSizes() const
-        {
-            return patchSizes_;
-        }
+            const wordList& patchNames() const
+            {
+                return patchNames_;
+            }
+
+            const labelList& patchSizes() const
+            {
+                return patchSizes_;
+            }
+
+            const labelList& patchStarts() const
+            {
+                return patchStarts_;
+            }
 
-        const labelList& patchStarts() const
-        {
-            return patchStarts_;
-        }
-        // Check
 
         // Edit
-        void createMesh();
-
-        // Write
 
+            //- Create the mesh
+            void createMesh();
 };
 
 

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

@@ -317,12 +317,12 @@ int main(int argc, char *argv[])
     }
 
     // Take over refinement levels and write to new time directory.
-    Pout<< "\nWriting extruded mesh to time = " << runTimeExtruded.timeName()
+    Info<< "\nWriting extruded mesh to time = " << runTimeExtruded.timeName()
         << nl << endl;
 
     mesh().write();
 
-    Pout<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }

applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict

@@ -69,8 +69,8 @@ structuredCoeffs
     // Renumber in columns (depthFirst) or in layers
     depthFirst true;
 
-    // Optional: reverse ordering
-    //reverse false;
+    // Reverse ordering
+    reverse false;
 }
 
 

applications/utilities/parallelProcessing/decomposePar/decomposeParDict

@@ -119,7 +119,10 @@ structuredCoeffs
 {
     // Patches to do 2D decomposition on. Structured mesh only; cells have
     // to be in 'columns' on top of patches.
-    patches     (bottomPatch);
+    patches     (movingWall);
+
+    // Method to use on the 2D subset
+    method      scotch;
 }
 
 //// Is the case distributed? Note: command-line argument -roots takes

applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H

@@ -2,31 +2,31 @@
 #define CGAL_PSURF_RINGS_H_
 
 // This file adapted from
-//     CGAL-3.7/examples/Jet_fitting_3//PolyhedralSurf_rings.h
-// Licensed under CGAL-3.7/LICENSE.FREE_USE
+//     CGAL-4.0/examples/Jet_fitting_3/PolyhedralSurf_rings.h
+// Licensed under CGAL-4.0/LICENSE.FREE_USE
 
 // Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
-// Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
-// Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
-// Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
-// Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
-// (Israel).  All rights reserved.
-
-// Permission is hereby granted, free of charge, to any person obtaining a copy
-// of this software and associated documentation files (the "Software"), to deal
-// in the Software without restriction, including without limitation the rights
-// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the Software is
-// furnished to do so, subject to the following conditions:
-
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
-// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-// SOFTWARE.
-
+// Utrecht University (The Netherlands),
+// ETH Zurich (Switzerland),
+// INRIA Sophia-Antipolis (France),
+// Max-Planck-Institute Saarbruecken (Germany),
+// and Tel-Aviv University (Israel).  All rights reserved.
+//
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the
+// "Software"), to deal in the Software without restriction, including
+// without limitation the rights to use, copy, modify, merge, publish,
+// distribute, sublicense, and/or sell copies of the Software, and to
+// permit persons to whom the Software is furnished to do so, subject to
+// the following conditions:
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
 #include <cassert>
 

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H

@@ -25,7 +25,8 @@ Class
     Foam::dynamicIndexedOctree
 
 Description
-    Non-pointer based hierarchical recursive searching
+    Non-pointer based hierarchical recursive searching.
+    Storage is dynamic, so elements can be deleted.
 
 SourceFiles
     dynamicIndexedOctree.C
@@ -451,6 +452,7 @@ public:
             );
         }
 
+
     // Member Functions
 
         // Access
@@ -656,6 +658,7 @@ public:
 
             label removeIndex(const label nodIndex, const label index);
 
+
         // Write
 
             //- Print tree. Either print all indices (printContent = true) or
@@ -671,6 +674,7 @@ public:
 
             void writeTreeInfo() const;
 
+
     // IOstream Operators
 
         friend Ostream& operator<< <Type>

src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,9 +27,6 @@ Class
 Description
     Simplified diagonal-based incomplete LU smoother for asymmetric matrices.
 
-    To improve efficiency, the residual is evaluated after every nSweeps
-    sweeps.
-
 SourceFiles
     DILUSmoother.C
 

src/OpenFOAM/meshes/meshShapes/face/face.C

@@ -318,6 +318,17 @@ int Foam::face::compare(const face& a, const face& b)
     {
         return 0;
     }
+    else if (sizeA == 1)
+    {
+        if (a[0] == b[0])
+        {
+            return 1;
+        }
+        else
+        {
+            return 0;
+        }
+    }
 
     const_circulator<face> aCirc(a);
     const_circulator<face> bCirc(b);
@@ -337,7 +348,7 @@ int Foam::face::compare(const face& a, const face& b)
     } while (bCirc.circulate(CirculatorBase::CLOCKWISE));
 
     // If the circulator has stopped then faces a and b do not share a matching
-    // point
+    // point. Doesn't work on matching, single element face.
     if (!bCirc.circulate())
     {
         return 0;

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -538,6 +538,9 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
     label cI
 )
 {
+    static label nWarnings = 0;
+    static const label maxWarnings = 100;
+
     const faceList& pFaces = mesh.faces();
     const labelList& pOwner = mesh.faceOwner();
 
@@ -553,19 +556,27 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
 
     if (tetBasePtI == -1)
     {
-        WarningIn
-        (
-            "Foam::List<Foam::tetIndices> "
-            "Foam::polyMeshTetDecomposition::faceTetIndices"
-            "("
-                "const polyMesh&, "
-                "label, "
-                "label"
-            ")"
-        )
-            << "No base point for face " << fI << ", " << f
-            << ", produces a valid tet decomposition."
-            << endl;
+        if (nWarnings < maxWarnings)
+        {
+            WarningIn
+            (
+                "Foam::List<Foam::tetIndices> "
+                "Foam::polyMeshTetDecomposition::faceTetIndices"
+                "("
+                    "const polyMesh&, "
+                    "label, "
+                    "label"
+                ")"
+            )   << "No base point for face " << fI << ", " << f
+                << ", produces a valid tet decomposition."
+                << endl;
+            nWarnings++;
+        }
+        if (nWarnings == maxWarnings)
+        {
+            Warning<< "Suppressing any further warnings." << endl;
+            nWarnings++;
+        }
 
         tetBasePtI = 0;
     }

src/dynamicMesh/createShellMesh/createShellMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -137,7 +137,7 @@ public:
             //            be in top patch
             //      < 0 : face in opposite orientation as patch face. face will
             //            be in bottom patch
-            //      = 0 : for all side and internal faces
+            //      = 0 : for all side faces
             const labelList& faceToFaceMap() const
             {
                 return faceToFaceMap_;

src/dynamicMesh/polyMeshFilter/polyMeshFilter.C

@@ -35,6 +35,12 @@ License
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+namespace Foam
+{
+defineTypeNameAndDebug(polyMeshFilter, 0);
+}
+
+
 Foam::autoPtr<Foam::fvMesh> Foam::polyMeshFilter::copyMesh(const fvMesh& mesh)
 {
     polyTopoChange originalMeshToNewMesh(mesh);

src/dynamicMesh/polyMeshFilter/polyMeshFilter.H

@@ -25,6 +25,10 @@ Class
     Foam::polyMeshFilter
 
 Description
+    Filter the edges and faces of a polyMesh whilst satisfying the given mesh
+    quality criteria.
+
+    Works on a copy of the mesh.
 
 SourceFiles
     polyMeshFilter.C
@@ -102,10 +106,10 @@ class polyMeshFilter
         //  faces
         const scalar faceReductionFactor_;
 
-        //-
+        //- Maximum number of times a deleted point can be associated with the
+        //  creation of a bad face it is forced to be kept.
         const label maxPointErrorCount_;
 
-
         //- The minimum edge length for each edge
         scalarField minEdgeLen_;
 
@@ -195,6 +199,10 @@ class polyMeshFilter
 
 public:
 
+    //- Runtime type information
+    ClassName("polyMeshFilter");
+
+
     // Constructors
 
         //- Construct from fvMesh
@@ -225,6 +233,7 @@ public:
             //- Filter edges only.
             label filterEdges(const label nOriginalBadFaces);
 
+            //- Filter all faces that are in the face zone indirectPatchFaces
             label filterIndirectPatchFaces();
 };
 

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C

@@ -35,6 +35,12 @@ License
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+namespace Foam
+{
+defineTypeNameAndDebug(edgeCollapser, 0);
+}
+
+
 Foam::label Foam::edgeCollapser::longestEdge
 (
     const face& f,

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H

@@ -61,7 +61,7 @@ class face;
 class edge;
 
 /*---------------------------------------------------------------------------*\
-                           Class edgeCollapser Declaration
+                        Class edgeCollapser Declaration
 \*---------------------------------------------------------------------------*/
 
 class edgeCollapser
@@ -84,12 +84,18 @@ private:
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        //- Controls collapse of a face to an edge
         const scalar guardFraction_;
 
+        //- Only collapse face to a point if high aspect ratio
         const scalar maxCollapseFaceToPointSideLengthCoeff_;
 
+        //- Allow a face to be collapsed to a point early, before the test
+        //  to collapse to an edge
         const Switch allowEarlyCollapseToPoint_;
 
+        //- Fraction of maxCollapseFaceToPointSideLengthCoeff_ to use when
+        //  allowEarlyCollapseToPoint_ is on
         const scalar allowEarlyCollapseCoeff_;
 
 
@@ -266,8 +272,8 @@ public:
                 const dictionary& meshQualityDict
             );
 
-            // Check mesh and mark points on faces in error
-            // Returns boolList with points in error set
+            //- Check mesh and mark points on faces in error
+            //  Returns boolList with points in error set
             static label checkMeshQuality
             (
                 const polyMesh& mesh,
@@ -300,7 +306,7 @@ public:
                 polyTopoChange& meshMod
             ) const;
 
-            // Mark (in collapseEdge) any edges to collapse
+            //- Mark (in collapseEdge) any edges to collapse
             label markSmallEdges
             (
                 const scalarField& minEdgeLen,
@@ -309,7 +315,7 @@ public:
                 Map<point>& collapsePointToLocation
             ) const;
 
-            // Mark (in collapseEdge) any edges to merge
+            //- Mark (in collapseEdge) any edges to merge
             label markMergeEdges
             (
                 const scalar maxCos,
@@ -332,6 +338,7 @@ public:
                 Map<point>& collapsePointToLocation
             ) const;
 
+            //- Marks edges in the faceZone indirectPatchFaces for collapse
             void markIndirectPatchFaces
             (
                 PackedBoolList& collapseEdge,

src/finiteVolume/finiteVolume/fvc/fvcReconstruct.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,6 +56,48 @@ reconstruct
 
     const fvMesh& mesh = ssf.mesh();
 
+    surfaceVectorField faceVols
+    (
+        mesh.Sf()/(mesh.magSf()*mesh.nonOrthDeltaCoeffs())
+    );
+
+    faceVols.internalField() *= (1.0 -  mesh.weights().internalField());
+    forAll(faceVols.boundaryField(), patchi)
+    {
+        if (faceVols.boundaryField()[patchi].coupled())
+        {
+            faceVols.boundaryField()[patchi] *=
+                (1.0 -  mesh.weights().boundaryField()[patchi]);
+        }
+    }
+
+    tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
+    (
+        new GeometricField<GradType, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                "volIntegrate("+ssf.name()+')',
+                ssf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            inv(surfaceSum(mesh.Sf()*faceVols))&surfaceSum(faceVols*ssf),
+            zeroGradientFvPatchField<GradType>::typeName
+        )
+    );
+
+    treconField().correctBoundaryConditions();
+
+    return treconField;
+}
+/*
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+
+    const fvMesh& mesh = ssf.mesh();
+
     tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
     (
         new GeometricField<GradType, fvPatchField, volMesh>
@@ -78,6 +120,7 @@ reconstruct
 
     return treconField;
 }
+*/
 
 
 template<class Type>

src/fvOptions/Make/files

@@ -28,11 +28,12 @@ $(derivedSources)/rotorDiskSource/trimModel/trimModel/trimModelNew.C
 $(derivedSources)/rotorDiskSource/trimModel/fixed/fixedTrim.C
 $(derivedSources)/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
 
-interRegion = $(derivedSources)/interRegionHeatTransferModel
-$(interRegion)/constantHeatTransfer/constantHeatTransfer.C
-$(interRegion)/interRegionHeatTransferModel/interRegionHeatTransferModel.C
-$(interRegion)/tabulatedHeatTransfer/tabulatedHeatTransfer.C
-$(interRegion)/variableHeatTransfer/variableHeatTransfer.C
+interRegion = sources/interRegion
+$(interRegion)/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
+$(interRegion)/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
+$(interRegion)/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
+$(interRegion)/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
+$(interRegion)/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
 
 
 /* constraints */

src/fvOptions/fvOptions/fvIOoptionList.C

@@ -45,14 +45,15 @@ Foam::IOobject Foam::fv::IOoptionList::createIOobject
 
     if (io.headerOk())
     {
-        Info<< "Creating field source list from " << io.name() << nl << endl;
+        Info<< "Creating fintite volume options from " << io.name() << nl
+            << endl;
 
         io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
         return io;
     }
     else
     {
-        Info<< "No field sources present" << nl << endl;
+        Info<< "No finite volume options present" << nl << endl;
 
         io.readOpt() = IOobject::NO_READ;
         return io;

src/fvOptions/fvOptions/fvOption.C

@@ -27,6 +27,7 @@ License
 #include "fvMesh.H"
 #include "fvMatrices.H"
 #include "volFields.H"
+#include "ListOps.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -86,9 +87,7 @@ void Foam::fv::option::setSelection(const dictionary& dict)
         }
         case smMapRegion:
         {
-            dict_.lookup("nbrModelName") >> nbrModelName_;
-            dict_.lookup("nbrRegionName") >> nbrRegionName_;
-            master_  = readBool(dict_.lookup("master"));
+            dict.lookup("nbrRegionName") >> nbrRegionName_;
             break;
         }
         case smAll:
@@ -131,7 +130,6 @@ void Foam::fv::option::setCellSet()
                     WarningIn("option::setCellSet()")
                         << "Unable to find owner cell for point " << points_[i]
                         << endl;
-
                 }
 
             }
@@ -167,40 +165,47 @@ void Foam::fv::option::setCellSet()
         }
         case smMapRegion:
         {
-            if (active_)
+            if (active_ && master_)
             {
                 IInfo<< "- selecting inter region mapping" << endl;
 
                 const fvMesh& nbrMesh =
                     mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
 
-                boundBox BB(mesh_.points(), false);
-                boundBox nbrBB(nbrMesh.points(), false);
-
-                if (nbrBB.overlaps(BB))
+                if (mesh_.name() == nbrMesh.name())
                 {
-                    // Dummy patches
-                    wordList cuttingPatches;
-                    HashTable<word> patchMap;
+                    FatalErrorIn("option::setCellIds()")
+                        << "Inter-region model selected, but local and "
+                        << "neighbour regions are the same: " << nl
+                        << "    local region: " << mesh_.name() << nl
+                        << "    secondary region: " << nbrMesh.name() << nl
+                        << exit(FatalError);
+                }
 
-                    secondaryToPrimaryInterpPtr_.reset
+                if (mesh_.bounds().overlaps(nbrMesh.bounds()))
+                {
+                    meshInterpPtr_.reset
                     (
-                        new meshToMesh
+                        new meshToMeshNew
                         (
-                            nbrMesh,
                             mesh_,
-                            patchMap,
-                            cuttingPatches
+                            nbrMesh,
+                            meshToMeshNew::interpolationMethodNames_.read
+                            (
+                                dict_.lookup("interpolationMethod")
+                            )
                         )
                     );
                 }
                 else
                 {
                     FatalErrorIn("option::setCellSet()")
-                        << "regions do not overlap " << nbrMesh.name()
-                        << " in region " << mesh_.name() << nl
+                        << "regions " << mesh_.name() << " and "
+                        << nbrMesh.name() <<  " do not intersect"
                         << exit(FatalError);
                 }
+
+                V_ = meshInterpPtr_->V();
             }
             break;
         }
@@ -244,7 +249,8 @@ Foam::fv::option::option
     const word& name,
     const word& modelType,
     const dictionary& dict,
-    const fvMesh& mesh
+    const fvMesh& mesh,
+    const bool master
 )
 :
     name_(name),
@@ -257,10 +263,9 @@ Foam::fv::option::option
     selectionMode_(selectionModeTypeNames_.read(dict_.lookup("selectionMode"))),
     cellSetName_("none"),
     V_(0.0),
-    secondaryToPrimaryInterpPtr_(),
-    nbrModelName_("none"),
+    meshInterpPtr_(),
     nbrRegionName_("none"),
-    master_(false),
+    master_(master),
     fieldNames_(),
     applied_()
 {
@@ -296,7 +301,7 @@ Foam::autoPtr<Foam::fv::option> Foam::fv::option::New
 {
     word modelType(coeffs.lookup("type"));
 
-    IInfo<< "Selecting source model type " << modelType << endl;
+    IInfo<< "Selecting finite volume options model type " << modelType << endl;
 
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(modelType);
@@ -317,12 +322,7 @@ Foam::autoPtr<Foam::fv::option> Foam::fv::option::New
 
 
 Foam::fv::option::~option()
-{
-    if (!secondaryToPrimaryInterpPtr_.empty())
-    {
-        secondaryToPrimaryInterpPtr_.clear();
-    }
-}
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -352,15 +352,7 @@ Foam::label Foam::fv::option::applyToField(const word& fieldName) const
         return 0;
     }
 
-    forAll(fieldNames_, i)
-    {
-        if (fieldNames_[i] == fieldName)
-        {
-            return i;
-        }
-    }
-
-    return -1;
+    return findIndex(fieldNames_, fieldName);
 }
 
 

src/fvOptions/fvOptions/fvOption.H

@@ -35,7 +35,7 @@ Description
         selectionMode   cellSet;                // cellSet // points //cellZone
                                                 // mapRegion
 Note:
-    On evaluation, source expects to be added to the rhs of the equation
+    On evaluation, source/sink options are to be added to the equation rhs
 
 SourceFiles
     fvOption.C
@@ -50,7 +50,7 @@ SourceFiles
 #include "volFieldsFwd.H"
 #include "cellSet.H"
 #include "autoPtr.H"
-#include "meshToMesh.H"
+#include "meshToMeshNew.H"
 
 #include "runTimeSelectionTables.H"
 
@@ -131,11 +131,8 @@ protected:
 
         // Data for smMapRegion only
 
-            //- Mesh to mesh mapping for map option
-            autoPtr<meshToMesh> secondaryToPrimaryInterpPtr_;
-
-            //- Name of the model in the neighbour mesh
-            word nbrModelName_;
+            //- Mesh to mesh interpolation object
+            autoPtr<meshToMeshNew> meshInterpPtr_;
 
             //- Name of the neighbour region to map
             word nbrRegionName_;
@@ -194,7 +191,8 @@ public:
             const word& name,
             const word& modelType,
             const dictionary& dict,
-            const fvMesh& mesh
+            const fvMesh& mesh,
+            const bool master = false
         );
 
         //- Return clone
@@ -288,15 +286,11 @@ public:
             //- Return const access to the total cell volume
             inline scalar V() const;
 
-            //- Return const access to the neighbour model name
-            inline const word& nbrModelName() const;
-
             //- Return const access to the neighbour region name
             inline const word& nbrRegionName() const;
 
-            //- Return const access to the mapToMap Ptr
-            inline const autoPtr<meshToMesh>
-                secondaryToPrimaryInterpPtr() const;
+            //- Return const access to the mapToMap pointer
+            inline const meshToMeshNew& meshInterp() const;
 
             //- Return const access to the cell set
             inline const labelList& cells() const;

src/fvOptions/fvOptions/fvOptionI.H

@@ -126,22 +126,22 @@ inline Foam::scalar& Foam::fv::option::duration()
 }
 
 
-inline const Foam::word& Foam::fv::option::nbrModelName() const
-{
-    return nbrModelName_;
-}
-
-
 inline const Foam::word& Foam::fv::option::nbrRegionName() const
 {
     return nbrRegionName_;
 }
 
 
-inline const Foam::autoPtr<Foam::meshToMesh>
-Foam::fv::option::secondaryToPrimaryInterpPtr() const
+inline const Foam::meshToMeshNew& Foam::fv::option::meshInterp() const
 {
-    return secondaryToPrimaryInterpPtr_;
+    if (!meshInterpPtr_.valid())
+    {
+        FatalErrorIn("const meshToMeshNew& meshInterp() const")
+            << "Interpolation object not set"
+            << abort(FatalError);
+    }
+
+    return meshInterpPtr_();
 }
 
 

src/fvOptions/fvOptions/fvOptionIO.C

@@ -90,8 +90,11 @@ void Foam::fv::option::writeData(Ostream& os) const
 bool Foam::fv::option::read(const dictionary& dict)
 {
     active_ = readBool(dict.lookup("active"));
-    timeStart_ = readScalar(dict.lookup("timeStart"));
-    duration_  = readScalar(dict.lookup("duration"));
+
+    if (dict.readIfPresent("timeStart", timeStart_))
+    {
+        dict.lookup("duration") >> duration_;
+    }
 
     coeffs_ = dict.subDict(type() + "Coeffs");
 

src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C

@@ -0,0 +1,230 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "interRegionExplicitPorositySource.H"
+#include "fvMesh.H"
+#include "fvMatrices.H"
+#include "porosityModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(interRegionExplicitPorositySource, 0);
+    addToRunTimeSelectionTable
+    (
+        option,
+        interRegionExplicitPorositySource,
+        dictionary
+    );
+}
+}
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::fv::interRegionExplicitPorositySource::initialise()
+{
+    if (!firstIter_)
+    {
+        return;
+    }
+
+    const word zoneName(name_ + ".porous");
+
+    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
+    const cellZoneMesh& cellZones = nbrMesh.cellZones();
+    label zoneID = cellZones.findZoneID(zoneName);
+
+    if (zoneID == -1)
+    {
+        cellZoneMesh& cz = const_cast<cellZoneMesh&>(cellZones);
+
+        zoneID = cz.size();
+
+        cz.setSize(zoneID + 1);
+
+        cz.set
+        (
+            zoneID,
+            new cellZone
+            (
+                zoneName,
+                nbrMesh.faceNeighbour(), // neighbour internal cells
+                zoneID,
+                cellZones
+            )
+        );
+
+        cz.clearAddressing();
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "void Foam::fv::interRegionExplicitPorositySource::initialise()"
+        )
+            << "Unable to create porous cellZone " << zoneName
+            << ": zone already exists"
+            << abort(FatalError);
+    }
+
+    porosityPtr_.reset
+    (
+        porosityModel::New
+        (
+            name_,
+            nbrMesh,
+            coeffs_,
+            zoneName
+        ).ptr()
+    ),
+
+    firstIter_ = false;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::interRegionExplicitPorositySource::interRegionExplicitPorositySource
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    option(name, modelType, dict, mesh, true),
+    porosityPtr_(NULL),
+    firstIter_(-1),
+    UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
+    rhoName_(coeffs_.lookupOrDefault<word>("rhoName", "rho")),
+    muName_(coeffs_.lookupOrDefault<word>("muName", "thermo:mu"))
+{
+    if (active_)
+    {
+        fieldNames_.setSize(1, UName_);
+        applied_.setSize(1, false);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::interRegionExplicitPorositySource::addSup
+(
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    initialise();
+
+    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
+
+    const volVectorField& U = eqn.psi();
+
+    volVectorField UNbr
+    (
+        IOobject
+        (
+            name_ + ".UNbr",
+            nbrMesh.time().timeName(),
+            nbrMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        nbrMesh,
+        dimensionedVector("zero", U.dimensions(), vector::zero)
+    );
+
+    // map local velocity onto neighbour region
+    meshInterp().mapSrcToTgt
+    (
+        U.internalField(),
+        plusEqOp<vector>(),
+        UNbr.internalField()
+    );
+
+    fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
+
+    if (eqn.dimensions() == dimForce)
+    {
+        const volScalarField& rhoNbr =
+            nbrMesh.lookupObject<volScalarField>(rhoName_);
+        const volScalarField& muNbr =
+            nbrMesh.lookupObject<volScalarField>(muName_);
+
+        porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr);
+    }
+    else
+    {
+        porosityPtr_->addResistance(nbrEqn);
+    }
+
+    // convert source from neighbour to local region
+    fvMatrix<vector> porosityEqn(U, eqn.dimensions());
+    scalarField& Udiag = porosityEqn.diag();
+    vectorField& Usource = porosityEqn.source();
+
+    Udiag.setSize(eqn.diag().size(), 0.0);
+    Usource.setSize(eqn.source().size(), vector::zero);
+
+    meshInterp().mapTgtToSrc(nbrEqn.diag(), plusEqOp<scalar>(), Udiag);
+    meshInterp().mapTgtToSrc(nbrEqn.source(), plusEqOp<vector>(), Usource);
+
+    eqn -= porosityEqn;
+}
+
+
+void Foam::fv::interRegionExplicitPorositySource::writeData(Ostream& os) const
+{
+    os  << indent << name_ << endl;
+    dict_.write(os);
+}
+
+
+bool Foam::fv::interRegionExplicitPorositySource::read(const dictionary& dict)
+{
+    if (option::read(dict))
+    {
+        coeffs_.readIfPresent("UName", UName_);
+        coeffs_.readIfPresent("rhoName", rhoName_);
+        coeffs_.readIfPresent("muName", muName_);
+
+        // reset the porosity model?
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::fv::interRegionExplicitPorositySource
+
+Description
+    Inter-region explicit porosity source
+
+    Sources described by, for example using the DarcyForchheimer model:
+
+        interRegionExplicitPorositySourceCoeffs
+        {
+            type            DarcyForchheimer;
+            DarcyForchheimerCoeffs
+            {
+                d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
+                f   f [0 -1 0 0 0 0 0] (0 0 0);
+
+                coordinateSystem
+                {
+                    e1  (0.70710678 0.70710678 0);
+                    e2  (0 0 1);
+                }
+            }
+        }
+
+Note:
+    The porous region must be selected as a cellZone.
+
+SourceFiles
+    interRegionExplicitPorositySource.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interRegionExplicitPorositySource_H
+#define interRegionExplicitPorositySource_H
+
+#include "fvOption.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class porosityModel;
+
+namespace fv
+{
+
+
+/*---------------------------------------------------------------------------*\
+                   Class interRegionExplicitPorositySource Declaration
+\*---------------------------------------------------------------------------*/
+
+class interRegionExplicitPorositySource
+:
+    public option
+{
+
+protected:
+
+    // Protected data
+
+        //- Run-time selectable porosity model
+        autoPtr<porosityModel> porosityPtr_;
+
+        //- First iteration
+        bool firstIter_;
+
+        //- Velocity field name, default = U
+        word UName_;
+
+        //- Density field name (compressible case only), default = rho
+        word rhoName_;
+
+        //- Dynamic viscosity field name (compressible case only)
+        //  default = thermo:mu
+        word muName_;
+
+
+    // Protected Member Functions
+
+        //- Initialise
+        void initialise();
+
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        interRegionExplicitPorositySource
+        (
+            const interRegionExplicitPorositySource&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const interRegionExplicitPorositySource&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("interRegionExplicitPorositySource");
+
+
+    // Constructors
+
+        //- Construct from components
+        interRegionExplicitPorositySource
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~interRegionExplicitPorositySource()
+    {}
+
+
+    // Member Functions
+
+        // Add explicit and implicit contributions
+
+            //- Vector
+            virtual void addSup
+            (
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
+
+
+        // I-O
+
+            //- Write data
+            virtual void writeData(Ostream&) const;
+
+            //- Read dictionary
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C

@@ -57,7 +57,7 @@ Foam::fv::constantHeatTransfer::constantHeatTransfer
     htcConst_(),
     AoV_()
 {
-    if (master_)
+    if (active() && master_)
     {
         htcConst_.reset
         (
@@ -91,19 +91,7 @@ Foam::fv::constantHeatTransfer::constantHeatTransfer
             )
         );
 
-        const DimensionedField<scalar, volMesh>& htcConsti =
-            htcConst_().dimensionedInternalField();
-        const DimensionedField<scalar, volMesh>& AoVi =
-            AoV_().dimensionedInternalField();
-        dimensionedScalar interVol
-        (
-            "V",
-            dimVolume,
-            secondaryToPrimaryInterpPtr_->V()
-        );
-
-        htc_.dimensionedInternalField() = htcConsti*AoVi*interVol/mesh.V();
-        htc_.correctBoundaryConditions();
+        htc_ = htcConst_()*AoV_();
     }
 }
 
@@ -116,10 +104,9 @@ Foam::fv::constantHeatTransfer::~constantHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::tmp<Foam::volScalarField>
-Foam::fv::constantHeatTransfer::calculateHtc()
+void Foam::fv::constantHeatTransfer::calculateHtc()
 {
-    return htc_;
+    // do nothing
 }
 
 

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H

@@ -87,7 +87,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual const tmp<volScalarField> calculateHtc();
+        virtual void calculateHtc();
 
 
         // I-O

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C

@@ -25,7 +25,7 @@ License
 
 #include "interRegionHeatTransferModel.H"
 #include "fluidThermo.H"
-#include "fvm.H"
+#include "fvmSup.H"
 #include "zeroGradientFvPatchFields.H"
 #include "fvcVolumeIntegrate.H"
 #include "fvOptionList.H"
@@ -41,10 +41,15 @@ namespace fv
 }
 
 
-// * * * * * * * * * * * *  Private member functions * * * * * * * * * * * //
+// * * * * * * * * * * * *  Protected member functions * * * * * * * * * * * //
 
-void Foam::fv::interRegionHeatTransferModel::check()
+void Foam::fv::interRegionHeatTransferModel::setNbrModel()
 {
+    if (!firstIter_)
+    {
+        return;
+    }
+
     const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
 
     const optionList& fvOptions = nbrMesh.lookupObject<optionList>("fvOptions");
@@ -66,11 +71,31 @@ void Foam::fv::interRegionHeatTransferModel::check()
 
     if (!nbrModelFound)
     {
-        FatalErrorIn("interRegionHeatTransferModel::check()")
-            << "Secondary source name not found" << nbrModelName_
+        FatalErrorIn("interRegionHeatTransferModel::setNbrModel()")
+            << "Neighbour model not found" << nbrModelName_
             << " in region " << nbrMesh.name() << nl
             << exit(FatalError);
     }
+
+    firstIter_ = false;
+}
+
+
+void Foam::fv::interRegionHeatTransferModel::correct()
+{
+    if (master_)
+    {
+        if (mesh_.time().timeIndex() != timeIndex_)
+        {
+            calculateHtc();
+            timeIndex_ = mesh_.time().timeIndex();
+        }
+    }
+    else
+    {
+        nbrModel().correct();
+        interpolate(nbrModel().htc(), htc_);
+    }
 }
 
 
@@ -84,9 +109,11 @@ Foam::fv::interRegionHeatTransferModel::interRegionHeatTransferModel
     const fvMesh& mesh
 )
 :
-    option(name, modelType, dict, mesh),
-    nbrModel_(),
+    option(name, modelType, dict, mesh, readBool(dict.lookup("master"))),
+    nbrModel_(NULL),
+    nbrModelName_(word::null),
     firstIter_(true),
+    timeIndex_(-1),
     htc_
     (
         IOobject
@@ -106,10 +133,14 @@ Foam::fv::interRegionHeatTransferModel::interRegionHeatTransferModel
         ),
         zeroGradientFvPatchScalarField::typeName
     ),
-    semiImplicit_(false)
+    semiImplicit_(false),
+    TName_(coeffs_.lookupOrDefault<word>("TName", "T")),
+    TNbrName_(coeffs_.lookupOrDefault<word>("TNbrName", "T"))
 {
     if (active())
     {
+        coeffs_.lookup("nbrModelName") >> nbrModelName_;
+
         coeffs_.lookup("fieldNames") >> fieldNames_;
         applied_.setSize(fieldNames_.size(), false);
 
@@ -132,19 +163,14 @@ void Foam::fv::interRegionHeatTransferModel::addSup
     const label fieldI
 )
 {
-    if (!secondaryToPrimaryInterpPtr_.valid())
-    {
-        return;
-    }
+    setNbrModel();
 
-    if (firstIter_)
-    {
-        check();
-        firstIter_ = false;
-    }
+    correct();
 
     const volScalarField& h = eqn.psi();
 
+    const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
+
     tmp<volScalarField> tTmapped
     (
         new volScalarField
@@ -157,8 +183,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            mesh_,
-            dimensionedScalar("T", dimTemperature, 0.0)
+            T
         )
     );
 
@@ -166,26 +191,10 @@ void Foam::fv::interRegionHeatTransferModel::addSup
 
     const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
 
-    const volScalarField& Tnbr = nbrMesh.lookupObject<volScalarField>("T");
+    const volScalarField& Tnbr =
+        nbrMesh.lookupObject<volScalarField>(TNbrName_);
 
-    secondaryToPrimaryInterpPtr_->interpolateInternalField
-    (
-        Tmapped,
-        Tnbr,
-        meshToMesh::MAP,
-        eqOp<scalar>()
-    );
-
-    if (!master_)
-    {
-        secondaryToPrimaryInterpPtr_->interpolateInternalField
-        (
-            htc_,
-            nbrModel_->calculateHtc(),
-            meshToMesh::CELL_VOLUME_WEIGHT,
-            eqOp<scalar>()
-        );
-    }
+    interpolate(Tnbr, Tmapped.internalField());
 
     if (debug)
     {
@@ -226,7 +235,6 @@ void Foam::fv::interRegionHeatTransferModel::addSup
     }
     else
     {
-        const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
         eqn += htc_*(Tmapped - T);
     }
 }

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H

@@ -26,7 +26,12 @@ Class
 
 Description
     Base class for inter region heat exchange. The derived classes must
-    provide the heat transfer coefficient (htc)
+    provide the heat transfer coeffisine (htc) which is used as follows
+    in the energy equation.
+
+     \f[
+        -htc*Tmapped + Sp(htc, T)
+     \f]
 
 \*---------------------------------------------------------------------------*/
 
@@ -60,27 +65,73 @@ private:
         //- Pointer to neighbour interRegionHeatTransferModel
         interRegionHeatTransferModel* nbrModel_;
 
+        //- Name of the model in the neighbour mesh
+        word nbrModelName_;
+
         //- First iteration
         bool firstIter_;
 
-
-    // Private members
-
-        //- Check coupled interRegionHeatTransferModel
-        void check();
+        //- Time index - used for updating htc
+        label timeIndex_;
 
 
 protected:
 
     // Protected data
 
-        //- Heat transfer coefficient [W/m2/k] by area/volume [1/m]
-        // registered on the master mesh
+        //- Heat transfer coefficient [W/m2/k] times area/volume [1/m]
         volScalarField htc_;
 
         //- Flag to activate semi-implicit coupling
         bool semiImplicit_;
 
+        //- Name of temperature field; default = "T"
+        word TName_;
+
+        //- Name of neighbour temperature field; default = "T"
+        word TNbrName_;
+
+
+    // Protected member functions
+
+        //- Set the neighbour interRegionHeatTransferModel
+        void setNbrModel();
+
+        //- Correct to calculate the inter-region heat transfer coefficient
+        void correct();
+
+        //- Interpolate field with nbrModel specified
+        template<class Type>
+        tmp<Field<Type> > interpolate
+        (
+            const interRegionHeatTransferModel& nbrModel,
+            const Field<Type>& field
+        ) const;
+
+        //- Interpolate field without nbrModel specified
+        template<class Type>
+        tmp<Field<Type> > interpolate
+        (
+            const Field<Type>& field
+        ) const;
+
+        //- Interpolate field with nbrModel specified
+        template<class Type>
+        void interpolate
+        (
+            const interRegionHeatTransferModel& nbrModel,
+            const Field<Type>& field,
+            Field<Type>& result
+        ) const;
+
+        //- Interpolate field without nbrModel specified
+        template<class Type>
+        void interpolate
+        (
+            const Field<Type>& field,
+            Field<Type>& result
+        ) const;
+
 
 public:
 
@@ -110,16 +161,19 @@ public:
         // Access
 
             //- Return the heat transfer coefficient
-            const volScalarField& htc() const
-            {
-                return  htc_;
-            }
+            inline const volScalarField& htc() const;
+
+            //- Return const access to the neighbour model
+            inline const interRegionHeatTransferModel& nbrModel() const;
+
+            //- Return access to the neighbour model
+            inline interRegionHeatTransferModel& nbrModel();
 
             //-Source term to fvMatrix<scalar>
             virtual void addSup(fvMatrix<scalar>& eqn, const label fieldI);
 
             //- Calculate heat transfer coefficient
-            virtual const tmp<volScalarField> calculateHtc() = 0;
+            virtual void calculateHtc() = 0;
 
 
         // I-O
@@ -138,6 +192,17 @@ public:
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "interRegionHeatTransferModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "interRegionHeatTransferModelTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelI.H

@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline const Foam::volScalarField&
+Foam::fv::interRegionHeatTransferModel::htc() const
+{
+    return  htc_;
+}
+
+
+inline const Foam::fv::interRegionHeatTransferModel&
+Foam::fv::interRegionHeatTransferModel::nbrModel() const
+{
+    if (nbrModel_ == NULL)
+    {
+        FatalErrorIn("const interRegionHeatTransferModel& nbrModel() const")
+            << "Neighbour model not set"
+            << abort(FatalError);
+    }
+
+    return *nbrModel_;
+}
+
+
+inline Foam::fv::interRegionHeatTransferModel&
+Foam::fv::interRegionHeatTransferModel::nbrModel()
+{
+    if (nbrModel_ == NULL)
+    {
+        FatalErrorIn("interRegionHeatTransferModel& nbrModel()")
+            << "Neighbour model not set"
+            << abort(FatalError);
+    }
+
+    return *nbrModel_;
+}
+
+
+// ************************************************************************* //

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelTemplates.C

@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::fv::interRegionHeatTransferModel::interpolate
+(
+    const interRegionHeatTransferModel& nbrModel,
+    const Field<Type>& field
+) const
+{
+    if (master_)
+    {
+        return meshInterp().mapTgtToSrc(field);
+    }
+    else
+    {
+        return (nbrModel.meshInterp().mapSrcToTgt(field));
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::fv::interRegionHeatTransferModel::interpolate
+(
+    const Field<Type>& field
+) const
+{
+    return interpolate(nbrModel(), field);
+}
+
+
+template<class Type>
+void Foam::fv::interRegionHeatTransferModel::interpolate
+(
+    const interRegionHeatTransferModel& nbrModel,
+    const Field<Type>& field,
+    Field<Type>& result
+) const
+{
+    if (master_)
+    {
+        meshInterp().mapTgtToSrc(field, plusEqOp<scalar>(), result);
+    }
+    else
+    {
+        nbrModel.meshInterp().mapSrcToTgt(field, plusEqOp<scalar>(), result);
+    }
+}
+
+
+template<class Type>
+void Foam::fv::interRegionHeatTransferModel::interpolate
+(
+    const Field<Type>& field,
+    Field<Type>& result
+) const
+{
+    return interpolate(nbrModel(), field, result);
+}
+
+
+// ************************************************************************* //

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C

@@ -94,6 +94,7 @@ Foam::fv::tabulatedHeatTransfer::tabulatedHeatTransfer
 :
     interRegionHeatTransferModel(name, modelType, dict, mesh),
     UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
+    UNbrName_(coeffs_.lookupOrDefault<word>("UNbrName", "U")),
     hTable_(),
     AoV_(),
     startTimeName_(mesh.time().timeName())
@@ -108,22 +109,16 @@ Foam::fv::tabulatedHeatTransfer::~tabulatedHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::tmp<Foam::volScalarField>
-Foam::fv::tabulatedHeatTransfer::calculateHtc()
+void Foam::fv::tabulatedHeatTransfer::calculateHtc()
 {
     const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName());
 
-    const volVectorField& UNbr = nbrMesh.lookupObject<volVectorField>(UName_);
+    const volVectorField& UNbr =
+        nbrMesh.lookupObject<volVectorField>(UNbrName_);
 
-    scalarField UMagNbrMapped(htc_.internalField().size(), 0.0);
+    const scalarField UMagNbr(mag(UNbr));
 
-    secondaryToPrimaryInterpPtr_->interpolateInternalField
-    (
-        UMagNbrMapped,
-        mag(UNbr),
-        meshToMesh::MAP,
-        eqOp<scalar>()
-    );
+    const scalarField UMagNbrMapped(interpolate(UMagNbr));
 
     const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
 
@@ -134,9 +129,7 @@ Foam::fv::tabulatedHeatTransfer::calculateHtc()
         htcc[i] = hTable()(mag(U[i]), UMagNbrMapped[i]);
     }
 
-    htcc = htcc*AoV()*secondaryToPrimaryInterpPtr_->V()/mesh_.V();
-
-    return htc_;
+    htcc = htcc*AoV();
 }
 
 

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H

@@ -60,6 +60,9 @@ private:
         //- Name of velocity field; default = U
         word UName_;
 
+        //- Name of neighbour velocity field; default = U
+        word UNbrName_;
+
         //- 2D look up table
         autoPtr<interpolation2DTable<scalar> > hTable_;
 
@@ -101,7 +104,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual const tmp<volScalarField> calculateHtc();
+        virtual void calculateHtc();
 
 
         // I-O

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C

@@ -55,7 +55,7 @@ Foam::fv::variableHeatTransfer::variableHeatTransfer
 )
 :
     interRegionHeatTransferModel(name, modelType, dict, mesh),
-    UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
+    UNbrName_(coeffs_.lookupOrDefault<word>("UNbrName", "U")),
     a_(0),
     b_(0),
     c_(0),
@@ -97,8 +97,7 @@ Foam::fv::variableHeatTransfer::~variableHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::tmp<Foam::volScalarField>
-Foam::fv::variableHeatTransfer::calculateHtc()
+void Foam::fv::variableHeatTransfer::calculateHtc()
 {
     const fvMesh& nbrMesh =
         mesh_.time().lookupObject<fvMesh>(nbrRegionName());
@@ -109,26 +108,18 @@ Foam::fv::variableHeatTransfer::calculateHtc()
     const fluidThermo& nbrThermo =
         nbrMesh.lookupObject<fluidThermo>("thermophysicalProperties");
 
-    const volVectorField& U = nbrMesh.lookupObject<volVectorField>(UName_);
+    const volVectorField& UNbr =
+        nbrMesh.lookupObject<volVectorField>(UNbrName_);
 
-    const volScalarField Re(mag(U)*ds_*nbrThermo.rho()/nbrTurb.mut());
+    const volScalarField ReNbr(mag(UNbr)*ds_*nbrThermo.rho()/nbrTurb.mut());
 
-    const volScalarField Nu(a_*pow(Re, b_)*pow(Pr_, c_));
+    const volScalarField NuNbr(a_*pow(ReNbr, b_)*pow(Pr_, c_));
 
-    scalarField htcNbrMapped(htc_.internalField().size(), 0.0);
+    const scalarField htcNbr(NuNbr*nbrTurb.kappaEff()/ds_);
 
-    secondaryToPrimaryInterpPtr_->interpolateInternalField
-    (
-        htcNbrMapped,
-        Nu*nbrTurb.kappaEff()/ds_,
-        meshToMesh::MAP,
-        eqOp<scalar>()
-    );
+    const scalarField htcNbrMapped(interpolate(htcNbr));
 
-    htc_.internalField() =
-        htcNbrMapped*AoV_*secondaryToPrimaryInterpPtr_->V()/mesh_.V();
-
-    return htc_;
+    htc_.internalField() = htcNbrMapped*AoV_;
 }
 
 
@@ -150,7 +141,7 @@ bool Foam::fv::variableHeatTransfer::read(const dictionary& dict)
 {
     if (option::read(dict))
     {
-        coeffs_.readIfPresent("UName", UName_);
+        coeffs_.readIfPresent("UNbrName", UNbrName_);
 
         coeffs_.readIfPresent("a", a_);
         coeffs_.readIfPresent("b", b_);

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H

@@ -65,8 +65,8 @@ private:
 
     // Private data
 
-        //- Name of velocity field; default = U
-        word UName_;
+        //- Name of neighbour velocity field; default = U
+        word UNbrName_;
 
         //- Model constants
         scalar a_;
@@ -108,7 +108,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual const tmp<volScalarField> calculateHtc();
+        virtual void calculateHtc();
 
 
         // I-O

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H

@@ -383,15 +383,25 @@ template<class CloudType>
 inline Foam::scalar Foam::ThermoCloud<CloudType>::Tmax() const
 {
     scalar T = -GREAT;
+    scalar n = 0;
     forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
         T = max(T, p.T());
+        n++;
     }
 
     reduce(T, maxOp<scalar>());
+    reduce(n, sumOp<label>());
 
-    return max(0.0, T);
+    if (n > 0)
+    {
+        return T;
+    }
+    else
+    {
+        return 0.0;
+    }
 }
 
 
@@ -399,15 +409,25 @@ template<class CloudType>
 inline Foam::scalar Foam::ThermoCloud<CloudType>::Tmin() const
 {
     scalar T = GREAT;
+    scalar n = 0;
     forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
         T = min(T, p.T());
+        n++;
     }
 
     reduce(T, minOp<scalar>());
+    reduce(n, sumOp<label>());
 
-    return max(0.0, T);
+    if (n > 0)
+    {
+        return T;
+    }
+    else
+    {
+        return 0.0;
+    }
 }
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -510,6 +510,9 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const CompositionModel<reactingCloudType>& composition =
         td.cloud().composition();
 
+    const scalar TMax = td.cloud().phaseChange().TMax(pc_, this->Tc_);
+    const scalar Tdash = min(T, TMax);
+    const scalar Tsdash = min(Ts, TMax);
 
     // Calculate mass transfer due to phase change
     td.cloud().phaseChange().calculate
@@ -520,8 +523,8 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         Pr,
         d,
         nus,
-        T,
-        Ts,
+        Tdash,
+        Tsdash,
         pc_,
         this->Tc_,
         YComponents,
@@ -541,7 +544,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         const label idc = composition.localToGlobalCarrierId(idPhase, i);
         const label idl = composition.globalIds(idPhase)[i];
 
-        const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
+        const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, Tdash);
         Sh -= dMassPC[i]*dh/dt;
     }
 
@@ -558,12 +561,12 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
             const label idc = composition.localToGlobalCarrierId(idPhase, i);
             const label idl = composition.globalIds(idPhase)[i];
 
-            const scalar Cp = composition.carrier().Cp(idc, pc_, Ts);
+            const scalar Cp = composition.carrier().Cp(idc, pc_, Tsdash);
             const scalar W = composition.carrier().W(idc);
             const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
 
             const scalar Dab =
-                composition.liquids().properties()[idl].D(pc_, Ts, Wc);
+                composition.liquids().properties()[idl].D(pc_, Tsdash, Wc);
 
             // Molar flux of species coming from the particle (kmol/m^2/s)
             N += Ni;

src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -404,20 +404,23 @@ void Foam::ParticleCollector<CloudType>::write()
     }
 
 
-    Field<scalar> faceMassTotal(mass_.size());
-    Field<scalar> faceMassFlowRate(massFlowRate_.size());
+    Field<scalar> faceMassTotal(mass_.size(), 0.0);
+    this->getModelProperty("massTotal", faceMassTotal);
+
+    Field<scalar> faceMassFlowRate(massFlowRate_.size(), 0.0);
+    this->getModelProperty("massFlowRate", faceMassFlowRate);
 
     forAll(faces_, faceI)
     {
         scalarList allProcMass(Pstream::nProcs());
         allProcMass[procI] = massTotal_[faceI];
         Pstream::gatherList(allProcMass);
-        faceMassTotal[faceI] = sum(allProcMass);
+        faceMassTotal[faceI] += sum(allProcMass);
 
         scalarList allProcMassFlowRate(Pstream::nProcs());
         allProcMassFlowRate[procI] = massFlowRate_[faceI];
         Pstream::gatherList(allProcMassFlowRate);
-        faceMassFlowRate[faceI] = sum(allProcMassFlowRate);
+        faceMassFlowRate[faceI] += sum(allProcMassFlowRate);
 
         Info<< "    face " << faceI
             << ": total mass = " << faceMassTotal[faceI]
@@ -470,20 +473,25 @@ void Foam::ParticleCollector<CloudType>::write()
 
     if (resetOnWrite_)
     {
-        forAll(faces_, faceI)
-        {
-            massFlowRate_[faceI] = 0.0;
-        }
+        Field<scalar> dummy(faceMassTotal.size(), 0.0);
+        this->setModelProperty("massTotal", dummy);
+        this->setModelProperty("massFlowRate", dummy);
+
         timeOld_ = timeNew;
         totalTime_ = 0.0;
     }
+    else
+    {
+        this->setModelProperty("massTotal", faceMassTotal);
+        this->setModelProperty("massFlowRate", faceMassFlowRate);
+    }
 
     forAll(faces_, faceI)
     {
         mass_[faceI] = 0.0;
+        massTotal_[faceI] = 0.0;
+        massFlowRate_[faceI] = 0.0;
     }
-
-    // writeProperties();
 }
 
 
@@ -552,8 +560,8 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
         (
             "Foam::ParticleCollector<CloudType>::ParticleCollector"
             "("
-                "const dictionary& dict,"
-                "CloudType& owner"
+                "const dictionary&,"
+                "CloudType&"
             ")"
         )
             << "Unknown mode " << mode << ".  Available options are "
@@ -565,8 +573,6 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     massFlowRate_.setSize(faces_.size(), 0.0);
 
     makeLogFile(faces_, points_, area_);
-
-    // readProperties(); AND initialise mass... fields
 }
 
 
@@ -579,6 +585,7 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     CloudFunctionObject<CloudType>(pc),
     mode_(pc.mode_),
     parcelType_(pc.parcelType_),
+    removeCollected_(pc.removeCollected_),
     points_(pc.points_),
     faces_(pc.faces_),
     faceTris_(pc.faceTris_),

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -200,8 +200,8 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
 (
     const label idc,
     const label idl,
-    const label p,
-    const label T
+    const scalar p,
+    const scalar T
 ) const
 {
     scalar dh = 0;
@@ -230,8 +230,8 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
                 "("
                     "const label, "
                     "const label, "
-                    "const label, "
-                    "const label"
+                    "const scalar, "
+                    "const scalar"
                 ") const"
             )   << "Unknown enthalpyTransfer type" << abort(FatalError);
         }
@@ -241,4 +241,21 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporation<CloudType>::TMax
+(
+    const scalar pIn,
+    const scalar TIn
+) const
+{
+    scalar T = -GREAT;
+    forAll(liquids_, i)
+    {
+        T = max(T, liquids_.properties()[i].pv(pIn, TIn));
+    }
+
+    return T;
+}
+
+
 // ************************************************************************* //

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,9 +127,12 @@ public:
         (
             const label idc,
             const label idl,
-            const label p,
-            const label T
+            const scalar p,
+            const scalar T
         ) const;
+
+        //- Return maximum/limiting temperature
+        virtual scalar TMax(const scalar pIn, const scalar TIn) const;
 };
 
 

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -293,8 +293,8 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
 (
     const label idc,
     const label idl,
-    const label p,
-    const label T
+    const scalar p,
+    const scalar T
 ) const
 {
     scalar dh = 0;
@@ -329,8 +329,8 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
                 "("
                     "const label, "
                     "const label, "
-                    "const label, "
-                    "const label"
+                    "const scalar, "
+                    "const scalar"
                 ") const"
             )   << "Unknown enthalpyTransfer type" << abort(FatalError);
         }
@@ -340,4 +340,21 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::TMax
+(
+    const scalar pIn,
+    const scalar TIn
+) const
+{
+    scalar T = -GREAT;
+    forAll(liquids_, i)
+    {
+        T = max(T, liquids_.properties()[i].pv(pIn, TIn));
+    }
+
+    return T;
+}
+
+
 // ************************************************************************* //

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -137,9 +137,12 @@ public:
         (
             const label idc,
             const label idl,
-            const label p,
-            const label T
+            const scalar p,
+            const scalar T
         ) const;
+
+        //- Return maximum/limiting temperature
+        virtual scalar TMax(const scalar pIn, const scalar TIn) const;
 };
 
 

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -170,14 +170,25 @@ Foam::scalar Foam::PhaseChangeModel<CloudType>::dh
 (
     const label idc,
     const label idl,
-    const label p,
-    const label T
+    const scalar p,
+    const scalar T
 ) const
 {
     return 0.0;
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::PhaseChangeModel<CloudType>::TMax
+(
+    const scalar,
+    const scalar
+) const
+{
+    return GREAT;
+}
+
+
 template<class CloudType>
 void Foam::PhaseChangeModel<CloudType>::addToPhaseChangeMass(const scalar dMass)
 {

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -180,10 +180,12 @@ public:
         (
             const label idc,
             const label idl,
-            const label p,
-            const label T
+            const scalar p,
+            const scalar T
         ) const;
 
+        //- Return maximum/limiting temperature
+        virtual scalar TMax(const scalar pIn, const scalar TIn) const;
 
         //- Add to phase change mass
         void addToPhaseChangeMass(const scalar dMass);

src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -572,7 +572,7 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
     regionModels::surfaceFilmModels::surfaceFilmModel& filmModel =
         const_cast<regionModels::surfaceFilmModels::surfaceFilmModel&>
         (
-            this->owner().db().objectRegistry::template
+            this->owner().db().time().objectRegistry::template
                 lookupObject<regionModels::surfaceFilmModels::surfaceFilmModel>
                 (
                     "surfaceFilmProperties"

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,9 +43,11 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const pointData& wDist)
     }
 }
 
+
 Foam::Istream& Foam::operator>>(Istream& is, pointData& wDist)
 {
     return is >> static_cast<pointEdgePoint&>(wDist) >> wDist.s_ >> wDist.v_;
 }
 
+
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,6 +148,13 @@ public:
                 TrackingData& td
             );
 
+    // Member Operators
+
+        // Needed for List IO
+        inline bool operator==(const pointData&) const;
+        inline bool operator!=(const pointData&) const;
+
+
     // IOstream Operators
 
         friend Ostream& operator<<(Ostream&, const pointData&);

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -215,4 +215,23 @@ inline bool Foam::pointData::updateEdge
 }
 
 
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::pointData::operator==(const Foam::pointData& rhs)
+const
+{
+    return
+        pointEdgePoint::operator==(rhs)
+     && (s() == rhs.s())
+     && (v() == rhs.v());
+}
+
+
+inline bool Foam::pointData::operator!=(const Foam::pointData& rhs)
+const
+{
+    return !(*this == rhs);
+}
+
+
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -258,6 +258,106 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 }
 
 
+void Foam::meshRefinement::testSyncPointList
+(
+    const string& msg,
+    const polyMesh& mesh,
+    const List<scalar>& fld
+)
+{
+    if (fld.size() != mesh.nPoints())
+    {
+        FatalErrorIn
+        (
+            "meshRefinement::testSyncPointList(const polyMesh&"
+            ", const List<scalar>&)"
+        )   << msg << endl
+            << "fld size:" << fld.size() << " mesh points:" << mesh.nPoints()
+            << abort(FatalError);
+    }
+
+    Pout<< "Checking field " << msg << endl;
+    scalarField minFld(fld);
+    syncTools::syncPointList
+    (
+        mesh,
+        minFld,
+        minEqOp<scalar>(),
+        GREAT
+    );
+    scalarField maxFld(fld);
+    syncTools::syncPointList
+    (
+        mesh,
+        maxFld,
+        maxEqOp<scalar>(),
+        -GREAT
+    );
+    forAll(minFld, pointI)
+    {
+        const scalar& minVal = minFld[pointI];
+        const scalar& maxVal = maxFld[pointI];
+        if (mag(minVal-maxVal) > SMALL)
+        {
+            Pout<< msg << " at:" << mesh.points()[pointI] << nl
+                << "    minFld:" << minVal << nl
+                << "    maxFld:" << maxVal << nl
+                << endl;
+        }
+    }
+}
+
+
+void Foam::meshRefinement::testSyncPointList
+(
+    const string& msg,
+    const polyMesh& mesh,
+    const List<point>& fld
+)
+{
+    if (fld.size() != mesh.nPoints())
+    {
+        FatalErrorIn
+        (
+            "meshRefinement::testSyncPointList(const polyMesh&"
+            ", const List<point>&)"
+        )   << msg << endl
+            << "fld size:" << fld.size() << " mesh points:" << mesh.nPoints()
+            << abort(FatalError);
+    }
+
+    Pout<< "Checking field " << msg << endl;
+    pointField minFld(fld);
+    syncTools::syncPointList
+    (
+        mesh,
+        minFld,
+        minMagSqrEqOp<point>(),
+        point(GREAT, GREAT, GREAT)
+    );
+    pointField maxFld(fld);
+    syncTools::syncPointList
+    (
+        mesh,
+        maxFld,
+        maxMagSqrEqOp<point>(),
+        vector::zero
+    );
+    forAll(minFld, pointI)
+    {
+        const point& minVal = minFld[pointI];
+        const point& maxVal = maxFld[pointI];
+        if (mag(minVal-maxVal) > SMALL)
+        {
+            Pout<< msg << " at:" << mesh.points()[pointI] << nl
+                << "    minFld:" << minVal << nl
+                << "    maxFld:" << maxVal << nl
+                << endl;
+        }
+    }
+}
+
+
 void Foam::meshRefinement::checkData()
 {
     Pout<< "meshRefinement::checkData() : Checking refinement structure."

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -879,6 +879,20 @@ public:
             //- Debugging: check that all faces still obey start()>end()
             void checkData();
 
+            static void testSyncPointList
+            (
+                const string& msg,
+                const polyMesh& mesh,
+                const List<scalar>& fld
+            );
+
+            static void testSyncPointList
+            (
+                const string& msg,
+                const polyMesh& mesh,
+                const List<point>& fld
+            );
+
             //- Compare two lists over all boundary faces
             template<class T>
             void testSyncBoundaryFaceList

src/meshTools/Make/files

@@ -201,4 +201,6 @@ regionCoupled/GAMG/interfaces/regionCoupledGAMGInterface/regionCoupledWallGAMGIn
 regionCoupled/GAMG/interfaceFields/regionCoupledGAMGInterfaceField/regionCoupledGAMGInterfaceField.C
 regionCoupled/GAMG/interfaceFields/regionCoupledGAMGInterfaceField/regionCoupledWallGAMGInterfaceField.C
 
+tetOverlapVolume/tetOverlapVolume.C
+
 LIB = $(FOAM_LIBBIN)/libmeshTools

src/meshTools/algorithms/PointEdgeWave/pointEdgePoint.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,8 +25,8 @@ Class
     Foam::pointEdgePoint
 
 Description
-    Holds information regarding nearest wall point. Used in pointEdgeWave.
-    (so not standard meshWave)
+    Holds information regarding nearest wall point. Used in PointEdgeWave.
+    (so not standard FaceCellWave)
     To be used in wall distance calculation.
 
 SourceFiles
@@ -116,7 +116,7 @@ public:
             inline scalar distSqr() const;
 
 
-        // Needed by meshWave
+        // Needed by PointEdgeWave
 
             //- Check whether origin has been changed at all or
             //  still contains original (invalid) value.

src/meshTools/algorithms/PointEdgeWave/pointEdgePointI.H

@@ -310,14 +310,14 @@ inline bool Foam::pointEdgePoint::equal
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::pointEdgePoint::operator==(const Foam::pointEdgePoint& rhs)
- const
+const
 {
-    return origin() == rhs.origin();
+    return (origin() == rhs.origin()) && (distSqr() == rhs.distSqr());
 }
 
 
 inline bool Foam::pointEdgePoint::operator!=(const Foam::pointEdgePoint& rhs)
- const
+const
 {
     return !(*this == rhs);
 }

src/meshTools/tetOverlapVolume/tetOverlapVolume.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,9 +54,9 @@ Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
     const tetPoints& tetB
 ) const
 {
-    tetPointRef::tetIntersectionList insideTets;
+    static tetPointRef::tetIntersectionList insideTets;
     label nInside = 0;
-    tetPointRef::tetIntersectionList cutInsideTets;
+    static tetPointRef::tetIntersectionList cutInsideTets;
     label nCutInside = 0;
 
     tetPointRef::storeOp inside(insideTets, nInside);
@@ -130,6 +130,135 @@ Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
 
 // * * * * * * * * * * * Public Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::tetOverlapVolume::cellCellOverlapMinDecomp
+(
+    const primitiveMesh& meshA,
+    const label cellAI,
+    const primitiveMesh& meshB,
+    const label cellBI,
+    const treeBoundBox& cellBbB,
+    const scalar threshold
+) const
+{
+    const cell& cFacesA = meshA.cells()[cellAI];
+    const point& ccA = meshA.cellCentres()[cellAI];
+
+    const cell& cFacesB = meshB.cells()[cellBI];
+    const point& ccB = meshB.cellCentres()[cellBI];
+
+    scalar vol = 0.0;
+
+    forAll(cFacesA, cFA)
+    {
+        label faceAI = cFacesA[cFA];
+
+        const face& fA = meshA.faces()[faceAI];
+        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
+        if (!pyrA.overlaps(cellBbB))
+        {
+            continue;
+        }
+
+        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
+
+        label tetBasePtAI = 0;
+
+        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
+
+        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
+        {
+            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
+            label otherFacePtAI = fA.fcIndex(facePtAI);
+
+            label pt0I = -1;
+            label pt1I = -1;
+
+            if (ownA)
+            {
+                pt0I = fA[facePtAI];
+                pt1I = fA[otherFacePtAI];
+            }
+            else
+            {
+                pt0I = fA[otherFacePtAI];
+                pt1I = fA[facePtAI];
+            }
+
+            const tetPoints tetA
+            (
+                ccA,
+                tetBasePtA,
+                meshA.points()[pt0I],
+                meshA.points()[pt1I]
+            );
+            const treeBoundBox tetABb(tetA.bounds());
+
+
+            // Loop over tets of cellB
+            forAll(cFacesB, cFB)
+            {
+                label faceBI = cFacesB[cFB];
+
+                const face& fB = meshB.faces()[faceBI];
+                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
+                if (!pyrB.overlaps(pyrA))
+                {
+                    continue;
+                }
+
+                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
+
+                label tetBasePtBI = 0;
+
+                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
+
+                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
+                {
+                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
+                    label otherFacePtBI = fB.fcIndex(facePtBI);
+
+                    label pt0I = -1;
+                    label pt1I = -1;
+
+                    if (ownB)
+                    {
+                        pt0I = fB[facePtBI];
+                        pt1I = fB[otherFacePtBI];
+                    }
+                    else
+                    {
+                        pt0I = fB[otherFacePtBI];
+                        pt1I = fB[facePtBI];
+                    }
+
+                    const tetPoints tetB
+                    (
+                        ccB,
+                        tetBasePtB,
+                        meshB.points()[pt0I],
+                        meshB.points()[pt1I]
+                    );
+
+                    if (!tetB.bounds().overlaps(tetABb))
+                    {
+                        continue;
+                    }
+
+                    vol += tetTetOverlapVol(tetA, tetB);
+
+                    if (vol > threshold)
+                    {
+                        return true;
+                    }
+                }
+            }
+        }
+    }
+
+    return false;
+}
+
+
 Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
 (
     const primitiveMesh& meshA,

src/meshTools/tetOverlapVolume/tetOverlapVolume.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,7 @@ Class
     Foam::tetOverlapVolume
 
 Description
-    Calculates overlap volume of two tets.
+    Calculates the overlap volume of two cells using tetrahedral decomposition
 
 SourceFiles
     tetOverlapVolume.C
@@ -48,14 +48,14 @@ class polyMesh;
 class tetPoints;
 
 /*---------------------------------------------------------------------------*\
-                         Class tetOverlapVolume Declaration
+                      Class tetOverlapVolume Declaration
 \*---------------------------------------------------------------------------*/
 
 class tetOverlapVolume
 {
     // Private member functions
 
-        //- Tet Overlap Vol
+        //- Tet overlap volume
         scalar tetTetOverlapVol
         (
             const tetPoints& tetA,
@@ -94,6 +94,16 @@ public:
             const label cellBI
         ) const;
 
+        //- Return true if olverlap volume is greater than threshold
+        bool cellCellOverlapMinDecomp
+        (
+            const primitiveMesh& meshA,
+            const label cellAI,
+            const primitiveMesh& meshB,
+            const label cellBI,
+            const treeBoundBox& cellBbB,
+            const scalar threshold = 0.0
+        ) const;
 
         //- Calculates the overlap volume
         scalar cellCellOverlapVolumeMinDecomp
@@ -112,7 +122,6 @@ public:
 
 } // End namespace Foam
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/parallel/decompose/decompositionMethods/structuredDecomp/topoDistanceDataI.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -124,8 +124,7 @@ inline bool Foam::topoDistanceData::updateCell
 {
     if (distance_ == -1)
     {
-        data_ = neighbourInfo.data_;
-        distance_ = neighbourInfo.distance_ + 1;
+        operator=(neighbourInfo);
         return true;
     }
     else
@@ -151,7 +150,8 @@ inline bool Foam::topoDistanceData::updateFace
 
     if (distance_ == -1)
     {
-        operator=(neighbourInfo);
+        data_ = neighbourInfo.data_;
+        distance_ = neighbourInfo.distance_ + 1;
         return true;
     }
     else

src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,6 +77,8 @@ Description
                 mean            on;
                 prime2Mean      on;
                 base            time;
+                window          10.0;
+                windowName      w1;
             }
             p
             {

src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,15 @@ void Foam::fieldAverage::addMeanField
 
         const word& fieldName = faItems_[fieldI].fieldName();
 
-        const word meanFieldName = fieldName + EXT_MEAN;
+        word meanFieldName = fieldName + EXT_MEAN;
+        if
+        (
+            (faItems_[fieldI].window() > 0)
+         && (faItems_[fieldI].windowName() != "")
+        )
+        {
+            meanFieldName = meanFieldName + "_" + faItems_[fieldI].windowName();
+        }
 
         Info<< "Reading/calculating field " << meanFieldName << nl << endl;
 
@@ -100,7 +108,16 @@ void Foam::fieldAverage::addPrime2MeanField
 
         const word& fieldName = faItems_[fieldI].fieldName();
 
-        const word meanFieldName = fieldName + EXT_PRIME2MEAN;
+        word meanFieldName = fieldName + EXT_PRIME2MEAN;
+        if
+        (
+            (faItems_[fieldI].window() > 0)
+         && (faItems_[fieldI].windowName() != "")
+        )
+        {
+            meanFieldName = meanFieldName + "_" + faItems_[fieldI].windowName();
+        }
+
         Info<< "Reading/calculating field " << meanFieldName << nl << endl;
 
         if (obr_.foundObject<fieldType2>(meanFieldName))

src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,7 +54,8 @@ Foam::fieldAverageItem::fieldAverageItem()
     mean_(0),
     prime2Mean_(0),
     base_(ITER),
-    window_(-1.0)
+    window_(-1.0),
+    windowName_("")
 {}
 
 
@@ -64,7 +65,8 @@ Foam::fieldAverageItem::fieldAverageItem(const fieldAverageItem& faItem)
     mean_(faItem.mean_),
     prime2Mean_(faItem.prime2Mean_),
     base_(faItem.base_),
-    window_(faItem.window_)
+    window_(faItem.window_),
+    windowName_(faItem.windowName_)
 {}
 
 
@@ -94,6 +96,7 @@ void Foam::fieldAverageItem::operator=(const fieldAverageItem& rhs)
     prime2Mean_ = rhs.prime2Mean_;
     base_ = rhs.base_;
     window_ = rhs.window_;
+    windowName_ = rhs.windowName_;
 }
 
 

src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ Description
         prime2Mean      on;
         base            time; // iteration
         window          200;  // optional averaging window
+        windowName      w1;   // optional window name (default = "")
     }
     \endverbatim
 
@@ -107,6 +108,9 @@ private:
         //- Averaging window - defaults to -1 for 'all iters/time'
         scalar window_;
 
+        //- Averaging window name - defaults to 'window'
+        word windowName_;
+
 
 public:
 
@@ -171,6 +175,11 @@ public:
                 return window_;
             }
 
+            const word& windowName() const
+            {
+                return windowName_;
+            }
+
 
     // Member Operators
 
@@ -190,7 +199,8 @@ public:
              && a.mean_ == b.mean_
              && a.prime2Mean_ == b.prime2Mean_
              && a.base_ == b.base_
-             && a.window_ == b.window_;
+             && a.window_ == b.window_
+             && a.windowName_ == b.windowName_;
         }
 
         friend bool operator!=

src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,6 +46,7 @@ Foam::fieldAverageItem::fieldAverageItem(Istream& is)
     entry.lookup("prime2Mean") >> prime2Mean_;
     base_ = baseTypeNames_[entry.lookup("base")];
     window_ = entry.lookupOrDefault<scalar>("window", -1.0);
+    windowName_ = entry.lookupOrDefault<word>("windowName", "");
 }
 
 
@@ -66,6 +67,7 @@ Foam::Istream& Foam::operator>>(Istream& is, fieldAverageItem& faItem)
     entry.lookup("prime2Mean") >> faItem.prime2Mean_;
     faItem.base_ = faItem.baseTypeNames_[entry.lookup("base")];
     faItem.window_ = entry.lookupOrDefault<scalar>("window", -1.0);
+    faItem.windowName_ = entry.lookupOrDefault<word>("windowName", "");
 
     return is;
 }
@@ -90,6 +92,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const fieldAverageItem& faItem)
     {
         os.writeKeyword("window") << faItem.window_
             << token::END_STATEMENT << nl;
+
+        if (faItem.windowName_ != "")
+        {
+            os.writeKeyword("windowName") << faItem.windowName_
+                << token::END_STATEMENT << nl;
+        }
     }
 
     os  << token::END_BLOCK << nl;

src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -198,7 +198,7 @@ Foam::labelList Foam::structuredRenumber::renumber
         }
         else
         {
-            label layerI = cellData[cellI].distance()-1;
+            label layerI = cellData[cellI].distance();
             if (depthFirst_)
             {
                 orderedToOld[nLayers*cellData[cellI].data()+layerI] = cellI;

src/sampling/Make/files

@@ -59,8 +59,8 @@ $(meshToMesh)/meshToMesh.C
 $(meshToMesh)/calculateMeshToMeshAddressing.C
 $(meshToMesh)/calculateMeshToMeshWeights.C
 
-tetOverlapVolume = meshToMeshInterpolation/tetOverlapVolume
-$(tetOverlapVolume)/tetOverlapVolume.C
-
+meshToMeshNew = meshToMeshInterpolation/meshToMeshNew
+$(meshToMeshNew)/meshToMeshNew.C
+$(meshToMeshNew)/meshToMeshNewParallelOps.C
 
 LIB = $(FOAM_LIBBIN)/libsampling