diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 9e02c8de2e..21a7fac1e7 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@@ -1,16 +1,11 @@
 Info<< "Creating reaction model\n" << endl;
 
-autoPtr<combustionModels::psiChemistryCombustionModel> reaction
+autoPtr<combustionModels::psiCombustionModel> reaction
 (
-    combustionModels::psiChemistryCombustionModel::New
-    (
-        mesh
-    )
+    combustionModels::psiCombustionModel::New(mesh)
 );
 
-psiChemistryModel& chemistry = reaction->pChemistry();
-
-psiReactionThermo& thermo = chemistry.thermo();
+psiReactionThermo& thermo = reaction->thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index 5344eea3be..f2517529e2 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -31,7 +31,7 @@ Description
 
 #include "fvCFD.H"
 #include "turbulenceModel.H"
-#include "psiChemistryCombustionModel.H"
+#include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
 
diff --git a/applications/solvers/combustion/rhoReactingFoam/createFields.H b/applications/solvers/combustion/rhoReactingFoam/createFields.H
index 5d71e0f075..f58e6e71bd 100644
--- a/applications/solvers/combustion/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/rhoReactingFoam/createFields.H
@@ -1,16 +1,11 @@
 Info<< "Creating reaction model\n" << endl;
 
-autoPtr<combustionModels::rhoChemistryCombustionModel> reaction
+autoPtr<combustionModels::rhoCombustionModel> reaction
 (
-    combustionModels::rhoChemistryCombustionModel::New
-    (
-        mesh
-    )
+    combustionModels::rhoCombustionModel::New(mesh)
 );
 
-rhoChemistryModel& chemistry = reaction->pChemistry();
-
-rhoReactionThermo& thermo = chemistry.thermo();
+rhoReactionThermo& thermo = reaction->thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
diff --git a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
index aa8f1a5501..f3bab92f2f 100644
--- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "rhoChemistryCombustionModel.H"
+#include "rhoCombustionModel.H"
 #include "turbulenceModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
index 4c368e7cdd..6f41f328a3 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,7 @@ Description
 #include "fvCFD.H"
 #include "turbulenceModel.H"
 #include "basicReactingMultiphaseCloud.H"
-#include "rhoChemistryCombustionModel.H"
+#include "rhoCombustionModel.H"
 #include "radiationModel.H"
 #include "porousZones.H"
 #include "IObasicSourceList.H"
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H
index 25ceef7aa6..b78aeea132 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H
@@ -1,16 +1,11 @@
     Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<combustionModels::rhoChemistryCombustionModel> combustion
+    autoPtr<combustionModels::rhoCombustionModel> combustion
     (
-        combustionModels::rhoChemistryCombustionModel::New
-        (
-            mesh
-        )
+        combustionModels::rhoCombustionModel::New(mesh)
     );
 
-    rhoChemistryModel& chemistry = combustion->pChemistry();
-
-    rhoReactionThermo& thermo = chemistry.thermo();
+    rhoReactionThermo& thermo = combustion->thermo();
 
     SLGThermo slgThermo(mesh, thermo);
 
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index d6c61074e1..7b9dff602e 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,7 @@ Description
 #include "turbulenceModel.H"
 #include "basicThermoCloud.H"
 #include "coalCloud.H"
-#include "psiChemistryCombustionModel.H"
+#include "psiCombustionModel.H"
 #include "IObasicSourceList.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
index 9a2b17f1af..22451ca108 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
@@ -1,16 +1,11 @@
     Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<combustionModels::psiChemistryCombustionModel> combustion
+    autoPtr<combustionModels::psiCombustionModel> combustion
     (
-        combustionModels::psiChemistryCombustionModel::New
-        (
-            mesh
-        )
+        combustionModels::psiCombustionModel::New(mesh)
     );
 
-    psiChemistryModel& chemistry = combustion->pChemistry();
-
-    psiReactionThermo& thermo = chemistry.thermo();
+    psiReactionThermo& thermo = combustion->thermo();
 
     SLGThermo slgThermo(mesh, thermo);
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
index 52ed150808..1dfb6aa284 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
@@ -1,16 +1,11 @@
     Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<combustionModels::rhoChemistryCombustionModel> combustion
+    autoPtr<combustionModels::rhoCombustionModel> combustion
     (
-        combustionModels::rhoChemistryCombustionModel::New
-        (
-            mesh
-        )
+        combustionModels::rhoCombustionModel::New(mesh)
     );
 
-    rhoChemistryModel& chemistry = combustion->pChemistry();
-
-    rhoReactionThermo& thermo = chemistry.thermo();
+    rhoReactionThermo& thermo = combustion->thermo();
 
     SLGThermo slgThermo(mesh, thermo);
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
index 8d97ec99ff..61533fcd6b 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Application
     porousExplicitSourceReactingParcelFoam
 
 Description
-    Transient PISO solver for compressible, laminar or turbulent flow with
+    Transient PIMPLE solver for compressible, laminar or turbulent flow with
     reacting multiphase Lagrangian parcels for porous media, including explicit
     sources for mass, momentum and energy
 
@@ -42,7 +42,7 @@ Description
 #include "fvCFD.H"
 #include "turbulenceModel.H"
 #include "basicReactingMultiphaseCloud.H"
-#include "rhoChemistryCombustionModel.H"
+#include "rhoCombustionModel.H"
 #include "radiationModel.H"
 #include "porousZones.H"
 #include "IObasicSourceList.H"
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
index 433e3d0182..6b0dd8655f 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
@@ -1,16 +1,11 @@
     Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<combustionModels::psiChemistryCombustionModel> combustion
+    autoPtr<combustionModels::psiCombustionModel> combustion
     (
-        combustionModels::psiChemistryCombustionModel::New
-        (
-            mesh
-        )
+        combustionModels::psiCombustionModel::New(mesh)
     );
 
-    psiChemistryModel& chemistry = combustion->pChemistry();
-
-    psiReactionThermo& thermo = chemistry.thermo();
+    psiReactionThermo& thermo = combustion->thermo();
 
     SLGThermo slgThermo(mesh, thermo);
 
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index 06ef5ddac3..9badba12b7 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Application
     reactingParcelFilmFoam
 
 Description
-    Transient PISO solver for compressible, laminar or turbulent flow with
+    Transient PIMPLE solver for compressible, laminar or turbulent flow with
     reacting Lagrangian parcels, and surface film modelling.
 
 \*---------------------------------------------------------------------------*/
@@ -34,7 +34,7 @@ Description
 #include "turbulenceModel.H"
 #include "basicReactingCloud.H"
 #include "surfaceFilmModel.H"
-#include "psiChemistryCombustionModel.H"
+#include "psiCombustionModel.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/Make/options
index d632e22f58..ba327a029f 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@@ -20,12 +21,14 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
+    -lsampling \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
@@ -45,4 +48,5 @@ EXE_LIBS = \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
+    -lfieldSources \
     -lcombustionModels
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
index eac49bfc33..2ddaa2c301 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
@@ -10,8 +10,10 @@
 
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
-        solve(UEqn == -fvc::grad(p));
+        solve(UEqn == -fvc::grad(p) + sources(rho, U));
         K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
index 78bbe14bbf..6cc645bde5 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_h)")
+        mesh.divScheme("div(phi,Yi_hs)")
     )
 );
 
@@ -30,9 +30,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               ==
                 parcels.SYi(i, Yi)
               + combustion->R(Yi)
+              + sources(rho, Yi)
             );
 
             YiEqn.relax();
+
+            sources.constrain(YiEqn);
+
             YiEqn.solve(mesh.solver("Yi"));
 
             Yi.max(0.0);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
index 655510006d..3d27c50bc6 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
@@ -1,16 +1,11 @@
     Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<combustionModels::psiChemistryCombustionModel> combustion
+    autoPtr<combustionModels::psiCombustionModel> combustion
     (
-        combustionModels::psiChemistryCombustionModel::New
-        (
-            mesh
-        )
+        combustionModels::psiCombustionModel::New(mesh)
     );
 
-    psiChemistryModel& chemistry = combustion->pChemistry();
-
-    psiReactionThermo& thermo = chemistry.thermo();
+    psiReactionThermo& thermo = combustion->thermo();
 
     SLGThermo slgThermo(mesh, thermo);
 
@@ -102,3 +97,6 @@
         mesh,
         dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );
+
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
index 9107df21c7..c7c1ec0957 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
@@ -1,5 +1,5 @@
 {
-    fvScalarMatrix hEqn
+    fvScalarMatrix hsEqn
     (
         fvm::ddt(rho, hs)
       + mvConvection->fvmDiv(phi, hs)
@@ -10,11 +10,14 @@
       + parcels.Sh(hs)
       + radiation->Sh(thermo)
       + combustion->Sh()
+      + sources(rho, hs)
     );
 
-    hEqn.relax();
+    hsEqn.relax();
 
-    hEqn.solve();
+    sources.constrain(hsEqn);
+
+    hsEqn.solve();
 
     thermo.correct();
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index c504fdf2ea..511347cdce 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -2,7 +2,7 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
+HbyA = rAU*(UEqn == sources(rho, U))().H();
 
 if (pimple.transonic())
 {
@@ -25,6 +25,7 @@ if (pimple.transonic())
           - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
+          + sources(psi, p, rho.name())
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -56,6 +57,7 @@ else
           - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
+          + sources(psi, p, rho.name())
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -72,6 +74,7 @@ else
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 dpdt = fvc::ddt(p);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index a46cb30f80..9b0ccc1963 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Application
     reactingParcelFoam
 
 Description
-    Transient PISO solver for compressible, laminar or turbulent flow with
+    Transient PIMPLE solver for compressible, laminar or turbulent flow with
     reacting Lagrangian parcels.
 
 \*---------------------------------------------------------------------------*/
@@ -33,10 +33,11 @@ Description
 #include "fvCFD.H"
 #include "turbulenceModel.H"
 #include "basicReactingCloud.H"
-#include "psiChemistryCombustionModel.H"
+#include "psiCombustionModel.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H
index 59e3edd7b3..b6293f2c1f 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,13 +30,18 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    solve
+    fvScalarMatrix rhoEqn
     (
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
         parcels.Srho(rho)
+      + sources(rho)
     );
+
+    sources.constrain(rhoEqn);
+
+    rhoEqn.solve();
 }
 
 // ************************************************************************* //
diff --git a/applications/solvers/lagrangian/sprayFoam/Make/options b/applications/solvers/lagrangian/sprayFoam/Make/options
index 9d6d7d5398..647a53c1b4 100644
--- a/applications/solvers/lagrangian/sprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/Make/options
@@ -2,6 +2,7 @@ EXE_INC = \
     -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@@ -21,12 +22,14 @@ EXE_INC = \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude
 
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
+    -lsampling \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
@@ -47,4 +50,5 @@ EXE_LIBS = \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
+    -lfieldSources \
     -lcombustionModels
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options
index 772d4099b0..f32b47478d 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options
@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
@@ -23,12 +24,14 @@ EXE_INC = \
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude
 
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
+    -lsampling \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
@@ -50,4 +53,5 @@ EXE_LIBS = \
     -lengine \
     -lregionModels \
     -lsurfaceFilmModels \
+    -lfieldSources \
     -lcombustionModels
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
index 5c4fa0adc1..1bdc5cf996 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
@@ -2,7 +2,7 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
+HbyA = rAU*(UEqn == sources(rho, U))().H();
 
 if (pimple.transonic())
 {
@@ -25,6 +25,7 @@ if (pimple.transonic())
           - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
+          + sources(psi, p, rho.name())
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -56,6 +57,7 @@ else
           - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
+          + sources(psi, p, rho.name())
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -72,6 +74,7 @@ else
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 dpdt = fvc::ddt(p);
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
index a4ae147f51..baf1e40dc7 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,10 +35,11 @@ Description
 #include "engineMesh.H"
 #include "turbulenceModel.H"
 #include "basicSprayCloud.H"
-#include "psiChemistryCombustionModel.H"
+#include "psiCombustionModel.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -105,7 +106,7 @@ int main(int argc, char *argv[])
 
         if (runTime.write())
         {
-            chemistry.dQ()().write();
+            combustion->dQ()().write();
         }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
index 50e24a70df..fce272d042 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,10 +33,11 @@ Description
 #include "fvCFD.H"
 #include "turbulenceModel.H"
 #include "basicSprayCloud.H"
-#include "psiChemistryCombustionModel.H"
+#include "psiCombustionModel.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -98,7 +99,7 @@ int main(int argc, char *argv[])
 
         if (runTime.write())
         {
-            chemistry.dQ()().write();
+            combustion->dQ()().write();
         }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
index 910732fa7a..bcdd163834 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
@@ -78,10 +78,8 @@ SourceFiles
 namespace Foam
 {
 
-typedef treeDataPrimitivePatch<face, List, const pointField, point>
-    treeDataBPatch;
-
 typedef PrimitivePatch<face, List, const pointField, point> bPatch;
+typedef treeDataPrimitivePatch<bPatch> treeDataBPatch;
 
 /*---------------------------------------------------------------------------*\
                   Class backgroundMeshDecomposition Declaration
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 074d5b8ee3..25405d2ff3 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -154,6 +154,7 @@ castellatedMeshControls
             //faceZone sphere;
             //cellZone sphere;
             //cellZoneInside inside;  //outside/insidePoint
+
             //- Optional specification of what to do with faceZone faces:
             //      internal : keep them as internal faces (default)
             //      baffle   : create baffles from them. This gives more
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 395ee77cbe..e64f6cb78f 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -1018,6 +1018,11 @@ int main(int argc, char *argv[])
 
     // Update proc maps
     if (cellProcAddressing.headerOk())
+    if
+    (
+        cellProcAddressing.headerOk()
+     && cellProcAddressing.size() == mesh.nCells()
+    )
     {
         Info<< "Renumbering processor cell decomposition map "
             << cellProcAddressing.name() << endl;
@@ -1028,6 +1033,11 @@ int main(int argc, char *argv[])
         );
     }
     if (faceProcAddressing.headerOk())
+    if
+    (
+        faceProcAddressing.headerOk()
+     && faceProcAddressing.size() == mesh.nFaces()
+    )
     {
         Info<< "Renumbering processor face decomposition map "
             << faceProcAddressing.name() << endl;
@@ -1054,6 +1064,11 @@ int main(int argc, char *argv[])
         }
     }
     if (pointProcAddressing.headerOk())
+    if
+    (
+        pointProcAddressing.headerOk()
+     && pointProcAddressing.size() == mesh.nPoints()
+    )
     {
         Info<< "Renumbering processor point decomposition map "
             << pointProcAddressing.name() << endl;
@@ -1173,21 +1188,41 @@ int main(int argc, char *argv[])
 
     mesh.write();
     if (cellProcAddressing.headerOk())
+    if
+    (
+        cellProcAddressing.headerOk()
+     && cellProcAddressing.size() == mesh.nCells()
+    )
     {
         cellProcAddressing.instance() = mesh.facesInstance();
         cellProcAddressing.write();
     }
     if (faceProcAddressing.headerOk())
+    if
+    (
+        faceProcAddressing.headerOk()
+     && faceProcAddressing.size() == mesh.nFaces()
+    )
     {
         faceProcAddressing.instance() = mesh.facesInstance();
         faceProcAddressing.write();
     }
     if (pointProcAddressing.headerOk())
+    if
+    (
+        pointProcAddressing.headerOk()
+     && pointProcAddressing.size() == mesh.nPoints()
+    )
     {
         pointProcAddressing.instance() = mesh.facesInstance();
         pointProcAddressing.write();
     }
     if (boundaryProcAddressing.headerOk())
+    if
+    (
+        boundaryProcAddressing.headerOk()
+     && boundaryProcAddressing.size() == mesh.boundaryMesh().size()
+    )
     {
         boundaryProcAddressing.instance() = mesh.facesInstance();
         boundaryProcAddressing.write();
diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
index 4c2180f21d..1f58d59dde 100644
--- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict
+++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
@@ -238,6 +238,14 @@ FoamFile
 //        cos     0.01;       // Tolerance (max cos of angle)
 //    }
 //
+//    // Walk on faces in faceSet, starting from face nearest given position
+//    source  regionToFace;
+//    sourceInfo
+//    {
+//        set         f0;
+//        nearPoint   (0.1 0.1 0.005);
+//    }
+//
 //
 //
 // pointSet
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
index 5a2e599a64..1e49db7505 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
@@ -811,7 +811,6 @@ void Foam::vtkPV3Foam::renderPatchNames(vtkRenderer* renderer, const bool show)
 }
 
 
-
 void Foam::vtkPV3Foam::PrintSelf(ostream& os, vtkIndent indent) const
 {
     os  << indent << "Number of nodes: "
@@ -821,9 +820,10 @@ void Foam::vtkPV3Foam::PrintSelf(ostream& os, vtkIndent indent) const
         << (meshPtr_ ? meshPtr_->nCells() : 0) << "\n";
 
     os  << indent << "Number of available time steps: "
-        << (dbPtr_.valid() ? dbPtr_().times().size() : 0) << endl;
+        << (dbPtr_.valid() ? dbPtr_().times().size() : 0) << "\n";
 
     os  << indent << "mesh region: " << meshRegion_ << "\n";
 }
 
+
 // ************************************************************************* //
diff --git a/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C b/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C
index 0bc6773518..d02da85df5 100644
--- a/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C
+++ b/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,9 +41,9 @@ int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
     #include "setRootCase.H"
-#   include "createTime.H"
+    #include "createTime.H"
     instantList timeDirs = timeSelector::select0(runTime, args);
-#   include "createMesh.H"
+    #include "createMesh.H"
 
     forAll(timeDirs, timeI)
     {
@@ -90,7 +90,7 @@ int main(int argc, char *argv[])
             mesh
         );
 
-#       include "createPhi.H"
+        #include "createPhi.H"
 
         singlePhaseTransportModel laminarTransport(U, phi);
 
@@ -124,8 +124,8 @@ int main(int argc, char *argv[])
 
                 Info<< "Patch " << patchi
                     << " named " << currPatch.name()
-                    << " y+ : min: " << min(Yp) << " max: " << max(Yp)
-                    << " average: " << average(Yp) << nl << endl;
+                    << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
+                    << " average: " << gAverage(Yp) << nl << endl;
             }
         }
 
diff --git a/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C b/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C
index c8de884de9..7af0a901af 100644
--- a/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C
+++ b/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C
@@ -86,8 +86,8 @@ void calcIncompressibleYPlus
 
             Info<< "Patch " << patchi
                 << " named " << nutPw.patch().name()
-                << " y+ : min: " << min(Yp) << " max: " << max(Yp)
-                << " average: " << average(Yp) << nl << endl;
+                << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
+                << " average: " << gAverage(Yp) << nl << endl;
         }
     }
 
@@ -166,8 +166,8 @@ void calcCompressibleYPlus
 
             Info<< "Patch " << patchi
                 << " named " << mutPw.patch().name()
-                << " y+ : min: " << min(Yp) << " max: " << max(Yp)
-                << " average: " << average(Yp) << nl << endl;
+                << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
+                << " average: " << gAverage(Yp) << nl << endl;
         }
     }
 
diff --git a/applications/utilities/preProcessing/setFields/setFields.C b/applications/utilities/preProcessing/setFields/setFields.C
index 7b690d46f2..41c09b6053 100644
--- a/applications/utilities/preProcessing/setFields/setFields.C
+++ b/applications/utilities/preProcessing/setFields/setFields.C
@@ -92,7 +92,15 @@ bool setCellFieldType
                 field.boundaryField()[patchi].patchInternalField();
         }
 
-        field.write();
+        if (!field.write())
+        {
+            FatalErrorIn
+            (
+                "void setCellFieldType"
+                "(const fvMesh& mesh, const labelList& selectedCells,"
+                "Istream& fieldValueStream)"
+            ) << "Failed writing field " << fieldName << endl;
+        }
     }
     else
     {
@@ -102,6 +110,9 @@ bool setCellFieldType
             "(const fvMesh& mesh, const labelList& selectedCells,"
             "Istream& fieldValueStream)"
         ) << "Field " << fieldName << " not found" << endl;
+
+        // Consume value
+        (void)pTraits<Type>(fieldValueStream);
     }
 
     return true;
@@ -260,7 +271,15 @@ bool setFaceFieldType
             }
         }
 
-        field.write();
+        if (!field.write())
+        {
+            FatalErrorIn
+            (
+                "void setFaceFieldType"
+                "(const fvMesh& mesh, const labelList& selectedFaces,"
+                "Istream& fieldValueStream)"
+            )   << "Failed writing field " << field.name() << exit(FatalError);
+        }
     }
     else
     {
@@ -270,6 +289,9 @@ bool setFaceFieldType
             "(const fvMesh& mesh, const labelList& selectedFaces,"
             "Istream& fieldValueStream)"
         ) << "Field " << fieldName << " not found" << endl;
+
+        // Consume value
+        (void)pTraits<Type>(fieldValueStream);
     }
 
     return true;
diff --git a/src/combustionModels/FSD/FSD.C b/src/combustionModels/FSD/FSD.C
index c7b646d8bb..26f58505a7 100644
--- a/src/combustionModels/FSD/FSD.C
+++ b/src/combustionModels/FSD/FSD.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,7 +36,8 @@ namespace combustionModels
 template<class CombThermoType, class ThermoType>
 FSD<CombThermoType, ThermoType>::FSD
 (
-    const word& modelType, const fvMesh& mesh
+    const word& modelType,
+    const fvMesh& mesh
 )
 :
     singleStepCombustion<CombThermoType, ThermoType>(modelType, mesh),
@@ -62,14 +63,8 @@ FSD<CombThermoType, ThermoType>::FSD
         this->mesh(),
         dimensionedScalar("zero", dimless, 0.0)
     ),
-    YFuelFuelStream_
-    (
-        dimensionedScalar("YFuelStream", dimless, 1.0)
-    ),
-    YO2OxiStream_
-    (
-        dimensionedScalar("YOxiStream", dimless, 0.23)
-    ),
+    YFuelFuelStream_(dimensionedScalar("YFuelStream", dimless, 1.0)),
+    YO2OxiStream_(dimensionedScalar("YOxiStream", dimless, 0.23)),
     Cv_(readScalar(this->coeffs().lookup("Cv"))),
     C_(5.0),
     ftMin_(0.0),
@@ -91,15 +86,15 @@ FSD<CombThermoType, ThermoType>::~FSD()
 template<class CombThermoType, class ThermoType>
 void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
 {
-    this->singleMixture_.fresCorrect();
+    this->singleMixturePtr_->fresCorrect();
 
-    const label fuelI = this->singleMixture_.fuelIndex();
+    const label fuelI = this->singleMixturePtr_->fuelIndex();
 
-    const volScalarField& YFuel = this->thermo_->composition().Y()[fuelI];
+    const volScalarField& YFuel = this->thermoPtr_->composition().Y()[fuelI];
 
-    const volScalarField& YO2 = this->thermo_->composition().Y("O2");
+    const volScalarField& YO2 = this->thermoPtr_->composition().Y("O2");
 
-    const dimensionedScalar s = this->singleMixture_.s();
+    const dimensionedScalar s = this->singleMixturePtr_->s();
 
     ft_ =
         (s*YFuel - (YO2 - YO2OxiStream_))/(s*YFuelFuelStream_ + YO2OxiStream_);
@@ -195,7 +190,7 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
 
     volScalarField  deltaF
     (
-        lesModel.delta()/dimensionedScalar("flame",dimLength, 1.5e-3)
+        lesModel.delta()/dimensionedScalar("flame", dimLength, 1.5e-3)
     );
 
     // Linear correlation between delta and flame thickness
@@ -205,33 +200,33 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
 
     forAll(ft_, cellI)
     {
-        if(ft_[cellI] > ftMin_ && ft_[cellI] < ftMax_)
+        if (ft_[cellI] > ftMin_ && ft_[cellI] < ftMax_)
         {
             scalar ftCell = ft_[cellI];
 
-            if(ftVar[cellI] > ftVarMin_) //sub-grid beta pdf of ft_
+            if (ftVar[cellI] > ftVarMin_) //sub-grid beta pdf of ft_
             {
                 scalar ftVarc = ftVar[cellI];
                 scalar a =
                     max(ftCell*(ftCell*(1.0 - ftCell)/ftVarc - 1.0), 0.0);
                 scalar b = max(a/ftCell - a, 0.0);
 
-                for(int i=1; i<ftDim_; i++)
+                for (int i=1; i<ftDim_; i++)
                 {
                     scalar ft = i*deltaFt;
                     pc[cellI] += pow(ft, a-1.0)*pow(1.0 - ft, b - 1.0)*deltaFt;
                 }
 
-                for(int i=1; i<ftDim_; i++)
+                for (int i=1; i<ftDim_; i++)
                 {
                     scalar ft = i*deltaFt;
                     omegaFuelBar[cellI] +=
                         omegaFuel[cellI]/omegaF[cellI]
                        *exp
-                       (
-                        -sqr(ft - ftStoich)
-                        /(2.0*sqr(0.01*omegaF[cellI]))
-                       )
+                        (
+                           -sqr(ft - ftStoich)
+                           /(2.0*sqr(0.01*omegaF[cellI]))
+                        )
                        *pow(ft, a - 1.0)
                        *pow(1.0 - ft, b - 1.0)
                        *deltaFt;
@@ -241,31 +236,25 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
             else
             {
                 omegaFuelBar[cellI] =
-                    (omegaFuel[cellI]/omegaF[cellI])
-                   *exp
-                   (
-                    -sqr(ftCell - ftStoich)/(2.0*sqr(0.01*omegaF[cellI]))
-                   );
+                   omegaFuel[cellI]/omegaF[cellI]
+                  *exp(-sqr(ftCell - ftStoich)/(2.0*sqr(0.01*omegaF[cellI])));
             }
-
         }
         else
         {
-                omegaFuelBar[cellI] = 0.0;
+            omegaFuelBar[cellI] = 0.0;
         }
     }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-//  Combustion progress variable (c).
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // Combustion progress variable, c
 
     List<label> productsIndex(2, label(-1));
     {
         label i = 0;
-        forAll (this->singleMixture_.specieProd(), specieI)
+        forAll(this->singleMixturePtr_->specieProd(), specieI)
         {
-            if (this->singleMixture_.specieProd()[specieI] < 0)
+            if (this->singleMixturePtr_->specieProd()[specieI] < 0)
             {
                 productsIndex[i] = specieI;
                 i++;
@@ -276,9 +265,9 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
 
     // Flamelet probability of the progress c based on IFC (reuse pc)
     scalar YprodTotal = 0;
-    forAll (productsIndex, j)
+    forAll(productsIndex, j)
     {
-        YprodTotal += this->singleMixture_.Yprod0()[productsIndex[j]];
+        YprodTotal += this->singleMixturePtr_->Yprod0()[productsIndex[j]];
     }
 
     forAll(ft_, cellI)
@@ -312,10 +301,10 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
 
     volScalarField& products = tproducts();
 
-    forAll (productsIndex, j)
+    forAll(productsIndex, j)
     {
         label specieI = productsIndex[j];
-        const volScalarField& Yp = this->thermo_->composition().Y()[specieI];
+        const volScalarField& Yp = this->thermoPtr_->composition().Y()[specieI];
         products += Yp;
     }
 
@@ -326,7 +315,7 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
 
     pc = min(C_*c, scalar(1));
 
-    const volScalarField fres(this->singleMixture_.fres(fuelI));
+    const volScalarField fres(this->singleMixturePtr_->fres(fuelI));
 
     this->wFuel_ == mgft*pc*omegaFuelBar;
 }
diff --git a/src/combustionModels/FSD/FSDs.C b/src/combustionModels/FSD/FSDs.C
index 19facf7ab5..4ee6d7bb36 100644
--- a/src/combustionModels/FSD/FSDs.C
+++ b/src/combustionModels/FSD/FSDs.C
@@ -27,6 +27,9 @@ License
 
 #include "thermoPhysicsTypes.H"
 #include "psiCombustionModel.H"
+#include "psiThermoCombustion.H"
+#include "rhoCombustionModel.H"
+#include "rhoThermoCombustion.H"
 #include "FSD.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -38,10 +41,34 @@ namespace combustionModels
     makeCombustionTypesThermo
     (
         FSD,
-        psiCombustionModel,
-        gasThermoPhysics
+        psiThermoCombustion,
+        gasThermoPhysics,
+        psiCombustionModel
     );
 
+    makeCombustionTypesThermo
+    (
+        FSD,
+        psiThermoCombustion,
+        constGasThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        FSD,
+        rhoThermoCombustion,
+        gasThermoPhysics,
+        rhoCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        FSD,
+        rhoThermoCombustion,
+        constGasThermoPhysics,
+        rhoCombustionModel
+    );
 }
 }
 
diff --git a/src/combustionModels/Make/files b/src/combustionModels/Make/files
index 325fa1dde1..20de91ecc6 100644
--- a/src/combustionModels/Make/files
+++ b/src/combustionModels/Make/files
@@ -1,19 +1,17 @@
 combustionModel/combustionModel.C
 
-psiCombustionModel/psiCombustionModel.C
-psiCombustionModel/psiCombustionModelNew.C
+psiCombustionModel/psiCombustionModel/psiCombustionModel.C
+psiCombustionModel/psiCombustionModel/psiCombustionModelNew.C
+psiCombustionModel/psiThermoCombustion/psiThermoCombustion.C
+psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.C
 
-rhoCombustionModel/rhoCombustionModel.C
-rhoCombustionModel/rhoCombustionModelNew.C
+rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.C
+rhoCombustionModel/rhoCombustionModel/rhoCombustionModelNew.C
+rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.C
+rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
 
 infinitelyFastChemistry/infinitelyFastChemistrys.C
 
-psiChemistryCombustionModel/psiChemistryCombustionModel.C
-psiChemistryCombustionModel/psiChemistryCombustionModelNew.C
-
-rhoChemistryCombustionModel/rhoChemistryCombustionModel.C
-rhoChemistryCombustionModel/rhoChemistryCombustionModelNew.C
-
 PaSR/PaSRs.C
 
 FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
diff --git a/src/combustionModels/PaSR/PaSR.C b/src/combustionModels/PaSR/PaSR.C
index 16720de35b..cd12d896e0 100644
--- a/src/combustionModels/PaSR/PaSR.C
+++ b/src/combustionModels/PaSR/PaSR.C
@@ -25,29 +25,27 @@ License
 #include "PaSR.H"
 #include "fvmSup.H"
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class CombThermoType>
-Foam::combustionModels::PaSR<CombThermoType>::PaSR
+template<class Type>
+Foam::combustionModels::PaSR<Type>::PaSR
 (
     const word& modelType,
     const fvMesh& mesh
 )
 :
-    CombThermoType(modelType, mesh),
-    Cmix_(this->coeffs().lookup("Cmix")),
+    Type(modelType, mesh),
+    Cmix_(readScalar(this->coeffs().lookup("Cmix"))),
     turbulentReaction_(this->coeffs().lookup("turbulentReaction")),
     kappa_
     (
         IOobject
         (
-            "kappa",
+            "PaSR::kappa",
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
-            IOobject::AUTO_WRITE,
-            false
+            IOobject::AUTO_WRITE
         ),
         mesh,
         dimensionedScalar("kappa", dimless, 0.0)
@@ -63,24 +61,22 @@ Foam::combustionModels::PaSR<CombThermoType>::PaSR
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class CombThermoType>
-Foam::combustionModels::PaSR<CombThermoType>::~PaSR()
+template<class Type>
+Foam::combustionModels::PaSR<Type>::~PaSR()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-
-template<class CombThermoType>
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::PaSR<CombThermoType>::tc() const
+template<class Type>
+Foam::tmp<Foam::volScalarField> Foam::combustionModels::PaSR<Type>::tc() const
 {
-    return this->pChemistry_->tc();
+    return this->chemistryPtr_->tc();
 }
 
 
-template<class CombThermoType>
-void Foam::combustionModels::PaSR<CombThermoType>::correct()
+template<class Type>
+void Foam::combustionModels::PaSR<Type>::correct()
 {
     if (this->active())
     {
@@ -89,11 +85,11 @@ void Foam::combustionModels::PaSR<CombThermoType>::correct()
 
         if (!useReactionRate_)
         {
-            this->pChemistry_->solve(t - dt, dt);
+            this->chemistryPtr_->solve(t - dt, dt);
         }
         else
         {
-            this->pChemistry_->calculate();
+            this->chemistryPtr_->calculate();
         }
 
         if (turbulentReaction_)
@@ -108,20 +104,12 @@ void Foam::combustionModels::PaSR<CombThermoType>::correct()
             tmp<volScalarField> ttc(tc());
             const volScalarField& tc = ttc();
 
-            const dimensionedScalar e0
-            (
-                "e0",
-                sqr(dimLength)/pow3(dimTime),
-                SMALL
-            );
-
             forAll(epsilon, i)
             {
                 if (epsilon[i] > 0)
                 {
                     scalar tk =
-                        Cmix_.value()
-                       *Foam::sqrt(muEff[i]/rho[i]/(epsilon[i] + e0.value()));
+                        Cmix_*Foam::sqrt(muEff[i]/rho[i]/(epsilon[i] + SMALL));
 
                     // Chalmers PaSR model
                     if (!useReactionRate_)
@@ -148,9 +136,9 @@ void Foam::combustionModels::PaSR<CombThermoType>::correct()
 }
 
 
-template<class CombThermoType>
+template<class Type>
 Foam::tmp<Foam::fvScalarMatrix>
-Foam::combustionModels::PaSR<CombThermoType>::R(const volScalarField& Y) const
+Foam::combustionModels::PaSR<Type>::R(const volScalarField& Y) const
 {
     tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
 
@@ -160,16 +148,16 @@ Foam::combustionModels::PaSR<CombThermoType>::R(const volScalarField& Y) const
     {
         const label specieI = this->thermo().composition().species()[Y.name()];
 
-        Su += kappa_*this->pChemistry_->RR(specieI);
+        Su += kappa_*this->chemistryPtr_->RR(specieI);
     }
 
     return tSu;
 }
 
 
-template<class CombThermoType>
+template<class Type>
 Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::PaSR<CombThermoType>::dQ() const
+Foam::combustionModels::PaSR<Type>::dQ() const
 {
     tmp<volScalarField> tdQ
     (
@@ -193,16 +181,16 @@ Foam::combustionModels::PaSR<CombThermoType>::dQ() const
     if (this->active())
     {
         volScalarField& dQ = tdQ();
-        dQ = kappa_*this->pChemistry_->dQ();
+        dQ = kappa_*this->chemistryPtr_->dQ();
     }
 
     return tdQ;
 }
 
 
-template<class CombThermoType>
+template<class Type>
 Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::PaSR<CombThermoType>::Sh() const
+Foam::combustionModels::PaSR<Type>::Sh() const
 {
     tmp<volScalarField> tSh
     (
@@ -226,17 +214,17 @@ Foam::combustionModels::PaSR<CombThermoType>::Sh() const
     if (this->active())
     {
         scalarField& Sh = tSh();
-        Sh = kappa_*this->pChemistry_->Sh();
+        Sh = kappa_*this->chemistryPtr_->Sh();
     }
 
     return tSh;
 }
 
 
-template<class CombThermoType>
-bool Foam::combustionModels::PaSR<CombThermoType>::read()
+template<class Type>
+bool Foam::combustionModels::PaSR<Type>::read()
 {
-    if (CombThermoType::read())
+    if (Type::read())
     {
         this->coeffs().lookup("Cmix") >> Cmix_;
         this->coeffs().lookup("turbulentReaction") >> turbulentReaction_;
diff --git a/src/combustionModels/PaSR/PaSR.H b/src/combustionModels/PaSR/PaSR.H
index 125c9182c6..4c4c3347e5 100644
--- a/src/combustionModels/PaSR/PaSR.H
+++ b/src/combustionModels/PaSR/PaSR.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,15 +47,15 @@ namespace combustionModels
                             Class PaSR Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class CombThermoType>
+template<class Type>
 class PaSR
 :
-    public CombThermoType
+    public Type
 {
     // Private data
 
         //- Mixing constant
-        dimensionedScalar Cmix_;
+        scalar Cmix_;
 
         //- Turbulent reaction switch
         Switch turbulentReaction_;
@@ -88,11 +88,7 @@ public:
     // Constructors
 
         //- Construct from components
-        PaSR
-        (
-            const word& modelType,
-            const fvMesh& mesh
-        );
+        PaSR(const word& modelType, const fvMesh& mesh);
 
 
     //- Destructor
@@ -115,6 +111,7 @@ public:
             //-  Return source for enthalpy equation [kg/m/s3]
             virtual tmp<volScalarField> Sh() const;
 
+
     // I-O
 
             //- Update properties from given dictionary
diff --git a/src/combustionModels/PaSR/PaSRs.C b/src/combustionModels/PaSR/PaSRs.C
index 435fffb094..24978b55c4 100644
--- a/src/combustionModels/PaSR/PaSRs.C
+++ b/src/combustionModels/PaSR/PaSRs.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,8 +25,8 @@ License
 
 #include "makeCombustionTypes.H"
 
-#include "psiChemistryCombustionModel.H"
-#include "rhoChemistryCombustionModel.H"
+#include "psiChemistryCombustion.H"
+#include "rhoChemistryCombustion.H"
 #include "PaSR.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -35,19 +35,10 @@ namespace Foam
 {
 namespace combustionModels
 {
-    makeCombustionTypes
-    (
-        PaSR,
-        psiChemistryCombustionModel
-    );
-
-    makeCombustionTypes
-    (
-        PaSR,
-        rhoChemistryCombustionModel
-    );
+    makeCombustionTypes(PaSR, psiChemistryCombustion, psiCombustionModel);
+    makeCombustionTypes(PaSR, rhoChemistryCombustion, rhoCombustionModel);
 }
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// ************************************************************************* //
diff --git a/src/combustionModels/combustionModel/makeCombustionTypes.H b/src/combustionModels/combustionModel/makeCombustionTypes.H
index da6bbc3e3e..31f2e84006 100644
--- a/src/combustionModels/combustionModel/makeCombustionTypes.H
+++ b/src/combustionModels/combustionModel/makeCombustionTypes.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,45 +30,44 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#define makeCombustionTypesThermo(CombModel, Comb, Thermo)                    \
-                                                                              \
-    typedef CombModel<Comb, Thermo> CombModel##Comb##Thermo;                  \
+#define makeCombustionTypesThermo(CombModel, CombType, Thermo, Table)         \
                                                                               \
+    typedef CombModel<CombType, Thermo> CombModel##CombType##Thermo;          \
                                                                               \
     defineTemplateTypeNameAndDebugWithName                                    \
     (                                                                         \
-        CombModel##Comb##Thermo,                                              \
-        #CombModel"<"#Comb","#Thermo">",                                      \
+        CombModel##CombType##Thermo,                                          \
+        #CombModel"<"#CombType","#Thermo">",                                  \
         0                                                                     \
     );                                                                        \
                                                                               \
-                                                                              \
     addToRunTimeSelectionTable                                                \
     (                                                                         \
-        Comb,                                                                 \
-        CombModel##Comb##Thermo,                                              \
+        Table,                                                                \
+        CombModel##CombType##Thermo,                                          \
         dictionary                                                            \
     );
 
-#define makeCombustionTypes(CombModel, CombThermoType)                        \
+
+#define makeCombustionTypes(CombModel, CombType, Table)                       \
                                                                               \
-    typedef CombModel<CombThermoType>                                         \
-        CombModel##CombThermoType;                                            \
+    typedef CombModel<CombType> CombModel##CombType;                          \
                                                                               \
     defineTemplateTypeNameAndDebugWithName                                    \
     (                                                                         \
-        CombModel##CombThermoType,                                            \
-        #CombModel"<"#CombThermoType">",                                      \
+        CombModel##CombType,                                                  \
+        #CombModel"<"#CombType">",                                            \
         0                                                                     \
     );                                                                        \
                                                                               \
     addToRunTimeSelectionTable                                                \
     (                                                                         \
-        CombThermoType,                                                       \
-        CombModel##CombThermoType,                                            \
+        Table,                                                                \
+        CombModel##CombType,                                                  \
         dictionary                                                            \
     );
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.C b/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.C
index ceee9ce28d..558ec8b274 100644
--- a/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.C
+++ b/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,21 +59,22 @@ void infinitelyFastChemistry<CombThermoType, ThermoType>::correct()
 
     if (this->active())
     {
-        this->singleMixture_.fresCorrect();
+        this->singleMixturePtr_->fresCorrect();
 
-        const label fuelI = this->singleMixture_.fuelIndex();
+        const label fuelI = this->singleMixturePtr_->fuelIndex();
 
-        const volScalarField& YFuel = this->thermo_->composition().Y()[fuelI];
+        const volScalarField& YFuel =
+            this->thermoPtr_->composition().Y()[fuelI];
 
-        const dimensionedScalar s = this->singleMixture_.s();
+        const dimensionedScalar s = this->singleMixturePtr_->s();
 
-        if (this->thermo_->composition().contains("O2"))
+        if (this->thermoPtr_->composition().contains("O2"))
         {
-            const volScalarField& YO2 = this->thermo_->composition().Y("O2");
+            const volScalarField& YO2 = this->thermoPtr_->composition().Y("O2");
 
             this->wFuel_ ==
-                 this->rho()/(this->mesh().time().deltaT()*C_)
-                *min(YFuel, YO2/s.value());
+                this->rho()/(this->mesh().time().deltaT()*C_)
+               *min(YFuel, YO2/s.value());
         }
     }
 }
diff --git a/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H b/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H
index c17f058e83..778670ce4c 100644
--- a/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H
+++ b/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ SourceFiles
 
 #ifndef infinitelyFastChemistry_H
 #define infinitelyFastChemistry_H
+
 #include "singleStepCombustion.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -51,7 +52,7 @@ namespace combustionModels
 template<class CombThermoType, class ThermoType>
 class infinitelyFastChemistry
 :
-    public singleStepCombustion <CombThermoType, ThermoType>
+    public singleStepCombustion<CombThermoType, ThermoType>
 {
     // Private data
 
@@ -77,10 +78,7 @@ public:
     // Constructors
 
         //- Construct from components
-        infinitelyFastChemistry
-        (
-            const word& modelType, const fvMesh& mesh
-        );
+        infinitelyFastChemistry(const word& modelType, const fvMesh& mesh);
 
 
     //- Destructor
@@ -95,7 +93,7 @@ public:
             virtual void correct();
 
 
-    // I-O
+        // I-O
 
             //- Update properties
             virtual bool read();
diff --git a/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistrys.C b/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistrys.C
index 0dc939fd2b..b8a48ed0ce 100644
--- a/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistrys.C
+++ b/src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistrys.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,10 +26,9 @@ License
 #include "makeCombustionTypes.H"
 
 #include "thermoPhysicsTypes.H"
-#include "psiCombustionModel.H"
-#include "rhoCombustionModel.H"
+#include "psiThermoCombustion.H"
+#include "rhoThermoCombustion.H"
 #include "infinitelyFastChemistry.H"
-#include "singleStepCombustion.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -37,26 +36,36 @@ namespace Foam
 {
 namespace combustionModels
 {
-
     makeCombustionTypesThermo
     (
         infinitelyFastChemistry,
-        psiCombustionModel,
-        gasThermoPhysics
+        psiThermoCombustion,
+        gasThermoPhysics,
+        psiCombustionModel
     );
 
     makeCombustionTypesThermo
     (
         infinitelyFastChemistry,
-        psiCombustionModel,
-        constGasThermoPhysics
+        psiThermoCombustion,
+        constGasThermoPhysics,
+        psiCombustionModel
     );
 
     makeCombustionTypesThermo
     (
         infinitelyFastChemistry,
-        rhoCombustionModel,
-        gasThermoPhysics
+        rhoThermoCombustion,
+        gasThermoPhysics,
+        rhoCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        infinitelyFastChemistry,
+        rhoThermoCombustion,
+        constGasThermoPhysics,
+        rhoCombustionModel
     );
 }
 }
diff --git a/src/combustionModels/noCombustion/noCombustions.C b/src/combustionModels/noCombustion/noCombustions.C
index dba5597f1e..b5a986b9ed 100644
--- a/src/combustionModels/noCombustion/noCombustions.C
+++ b/src/combustionModels/noCombustion/noCombustions.C
@@ -27,6 +27,9 @@ License
 
 #include "psiCombustionModel.H"
 #include "rhoCombustionModel.H"
+#include "psiThermoCombustion.H"
+#include "rhoThermoCombustion.H"
+
 #include "noCombustion.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -38,12 +41,14 @@ namespace combustionModels
     makeCombustionTypes
     (
         noCombustion,
+        psiThermoCombustion,
         psiCombustionModel
     );
 
     makeCombustionTypes
     (
         noCombustion,
+        rhoThermoCombustion,
         rhoCombustionModel
     );
 }
diff --git a/src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModelNew.C b/src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModelNew.C
deleted file mode 100644
index 1ccba3fc8a..0000000000
--- a/src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModelNew.C
+++ /dev/null
@@ -1,78 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "psiChemistryCombustionModel.H"
-
-// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::combustionModels::psiChemistryCombustionModel>
-Foam::combustionModels::psiChemistryCombustionModel::New
-(
-    const fvMesh& mesh
-)
-{
-    const word combModelName
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "combustionProperties",
-                mesh.time().constant(),
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        ).lookup("combustionModel")
-    );
-
-    Info<< "Selecting combustion model " << combModelName << endl;
-
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(combModelName);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
-    {
-        FatalErrorIn
-        (
-            "psiChemistryCombustionModel::New"
-        )   << "Unknown psiChemistryCombustionModel type "
-            << combModelName << endl << endl
-            << "Valid  combustionModels are : " << endl
-            << dictionaryConstructorTablePtr_->toc()
-            << exit(FatalError);
-    }
-
-    const label tempOpen = combModelName.find('<');
-
-    const word className = combModelName(0, tempOpen);
-
-    return autoPtr<psiChemistryCombustionModel>
-        (cstrIter()(className, mesh));
-}
-
-
-// ************************************************************************* //
diff --git a/src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModel.C b/src/combustionModels/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.C
similarity index 60%
rename from src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModel.C
rename to src/combustionModels/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.C
index e78febd054..dc59186040 100644
--- a/src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModel.C
+++ b/src/combustionModels/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,54 +23,48 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "psiChemistryCombustionModel.H"
-
-/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
-
-namespace Foam
-{
-namespace combustionModels
-{
-    defineTypeNameAndDebug(psiChemistryCombustionModel, 0);
-    defineRunTimeSelectionTable(psiChemistryCombustionModel, dictionary);
-
-
+#include "psiChemistryCombustion.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-psiChemistryCombustionModel::psiChemistryCombustionModel
+Foam::combustionModels::psiChemistryCombustion::psiChemistryCombustion
 (
     const word& modelType,
     const fvMesh& mesh
 )
 :
-    combustionModel(modelType, mesh),
-    pChemistry_(psiChemistryModel::New(mesh))
+    psiCombustionModel(modelType, mesh),
+    chemistryPtr_(psiChemistryModel::New(mesh))
 {}
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-psiChemistryCombustionModel::~psiChemistryCombustionModel()
+Foam::combustionModels::psiChemistryCombustion::~psiChemistryCombustion()
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-
-bool psiChemistryCombustionModel::read()
+Foam::psiReactionThermo&
+Foam::combustionModels::psiChemistryCombustion::thermo()
 {
-    if (combustionModel::read())
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
+    return chemistryPtr_->thermo();
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-} // End namespace combustionModels
+const Foam::psiReactionThermo&
+Foam::combustionModels::psiChemistryCombustion::thermo() const
+{
+    return chemistryPtr_->thermo();
+}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::psiChemistryCombustion::rho() const
+{
+    return chemistryPtr_->thermo().rho();
+}
+
+
+// ************************************************************************* //
diff --git a/src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModel.H b/src/combustionModels/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.H
similarity index 51%
rename from src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModel.H
rename to src/combustionModels/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.H
index bde6be81b8..81c7a91139 100644
--- a/src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModel.H
+++ b/src/combustionModels/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,24 +22,21 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::psiChemistryCombustionModel
+    Foam::psiChemistryCombustion
 
 Description
-    Combustion models for compressibility-based thermodynamics
+    Compressibility-based chemistry model wrapper for combustion models
 
 SourceFiles
-    psiChemistryCombustionModelI.H
-    psiChemistryCombustionModel.C
-    psiChemistryCombustionModelNew.C
+    psiChemistryCombustion.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef psiChemistryCombustionModel_H
-#define psiChemistryCombustionModel_H
+#ifndef psiChemistryCombustion_H
+#define psiChemistryCombustion_H
 
-#include "combustionModel.H"
 #include "autoPtr.H"
-#include "runTimeSelectionTables.H"
+#include "psiCombustionModel.H"
 #include "psiChemistryModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -50,108 +47,59 @@ namespace combustionModels
 {
 
 /*---------------------------------------------------------------------------*\
-                     class psiChemistryCombustionModel Declaration
+                    class psiChemistryCombustion Declaration
 \*---------------------------------------------------------------------------*/
 
-class psiChemistryCombustionModel
+class psiChemistryCombustion
 :
-    public combustionModel
+    public psiCombustionModel
 {
-
     // Private Member Functions
 
         //- Construct as copy (not implemented)
-        psiChemistryCombustionModel(const psiChemistryCombustionModel&);
+        psiChemistryCombustion(const psiChemistryCombustion&);
 
         //- Disallow default bitwise assignment
-        void operator=(const psiChemistryCombustionModel&);
+        void operator=(const psiChemistryCombustion&);
 
 
 protected:
 
     // Protected data
 
-        //- Auto pointer to psiChemistry
-        autoPtr<psiChemistryModel> pChemistry_;
+        //- Pointer to chemistry model
+        autoPtr<psiChemistryModel> chemistryPtr_;
 
 
 public:
 
-    //- Runtime type information
-    TypeName("psiChemistryCombustionModel");
-
-
-    //- Declare run-time constructor selection tables
-    declareRunTimeSelectionTable
-    (
-        autoPtr,
-        psiChemistryCombustionModel,
-        dictionary,
-        (
-            const word& modelType,
-            const fvMesh& mesh
-        ),
-        (modelType, mesh)
-    );
-
-
     // Constructors
 
-
         //- Construct from components and thermo
-        psiChemistryCombustionModel
-        (
-            const word& modelType,
-            const fvMesh& mesh
-        );
-
-
-
-    //- Selector
-    static autoPtr<psiChemistryCombustionModel> New
-    (
-        const fvMesh& mesh
-    );
+        psiChemistryCombustion(const word& modelType, const fvMesh& mesh);
 
 
     //- Destructor
-    virtual ~psiChemistryCombustionModel();
+    virtual ~psiChemistryCombustion();
 
 
     // Member Functions
 
-
         //- Return access to the thermo package
-        inline psiChemistryModel& pChemistry();
+        virtual psiReactionThermo& thermo();
 
         //- Return const access to the thermo package
-        inline const psiChemistryModel& pChemistry() const;
+        virtual const psiReactionThermo& thermo() const;
 
-        //- Return const access to rho
-        inline tmp<volScalarField> rho() const;
-
-        //- Return const access to rho
-        inline const psiReactionThermo& thermo() const;
-
-        //- Return non const access to rho
-        inline psiReactionThermo& thermo();
-
-
-    // I-O
-
-        //- Update properties from given dictionary
-        virtual bool read();
+        //- Return const access to the density field
+        virtual tmp<volScalarField> rho() const;
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
 } // End namespace combustionModels
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "psiChemistryCombustionModelI.H"
+} // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/combustionModels/psiCombustionModel/psiCombustionModel.C b/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModel.C
similarity index 91%
rename from src/combustionModels/psiCombustionModel/psiCombustionModel.C
rename to src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModel.C
index 906934ab70..287f102101 100644
--- a/src/combustionModels/psiCombustionModel/psiCombustionModel.C
+++ b/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModel.C
@@ -38,24 +38,23 @@ namespace combustionModels
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-
 Foam::combustionModels::psiCombustionModel::psiCombustionModel
 (
     const word& modelType,
     const fvMesh& mesh
 )
 :
-    combustionModel(modelType, mesh),
-    thermo_(psiReactionThermo::New(mesh))
+    combustionModel(modelType, mesh)
 {}
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::combustionModels::psiCombustionModel::~psiCombustionModel()
 {}
 
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 bool Foam::combustionModels::psiCombustionModel::read()
 {
@@ -69,10 +68,5 @@ bool Foam::combustionModels::psiCombustionModel::read()
     }
 }
 
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::psiCombustionModel::rho() const
-{
-    return thermo_->rho();
-}
 
 // ************************************************************************* //
diff --git a/src/combustionModels/psiCombustionModel/psiCombustionModel.H b/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModel.H
similarity index 85%
rename from src/combustionModels/psiCombustionModel/psiCombustionModel.H
rename to src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModel.H
index e89e90d519..0206a5e70c 100644
--- a/src/combustionModels/psiCombustionModel/psiCombustionModel.H
+++ b/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModel.H
@@ -67,17 +67,6 @@ class psiCombustionModel
         void operator=(const psiCombustionModel&);
 
 
-protected:
-
-    // Protected data
-
-
-        //- Thermo package
-        autoPtr<psiReactionThermo> thermo_;
-
-
-
-
 public:
 
     //- Runtime type information
@@ -100,21 +89,13 @@ public:
 
     // Constructors
 
-
         //- Construct from components and thermo
-        psiCombustionModel
-        (
-            const word& modelType,
-            const fvMesh& mesh
-        );
+        psiCombustionModel(const word& modelType, const fvMesh& mesh);
 
 
 
     //- Selector
-    static autoPtr<psiCombustionModel> New
-    (
-        const fvMesh& mesh
-    );
+    static autoPtr<psiCombustionModel> New(const fvMesh& mesh);
 
 
     //- Destructor
@@ -124,13 +105,14 @@ public:
     // Member Functions
 
         //- Return access to the thermo package
-        inline psiReactionThermo& thermo();
+        virtual psiReactionThermo& thermo() = 0;
 
         //- Return const access to the thermo package
-        inline const psiReactionThermo& thermo() const;
+        virtual const psiReactionThermo& thermo() const = 0;
 
         //- Return tmp of rho
-        virtual tmp<volScalarField> rho() const;
+        virtual tmp<volScalarField> rho() const = 0;
+
 
     // I-O
 
@@ -141,12 +123,8 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
 } // End namespace combustionModels
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "psiCombustionModelI.H"
+} // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/combustionModels/psiCombustionModel/psiCombustionModelNew.C b/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModelNew.C
similarity index 100%
rename from src/combustionModels/psiCombustionModel/psiCombustionModelNew.C
rename to src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModelNew.C
diff --git a/src/combustionModels/psiCombustionModel/psiCombustionModels.C b/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C
similarity index 100%
rename from src/combustionModels/psiCombustionModel/psiCombustionModels.C
rename to src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C
diff --git a/src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModel.C b/src/combustionModels/psiCombustionModel/psiThermoCombustion/psiThermoCombustion.C
similarity index 60%
rename from src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModel.C
rename to src/combustionModels/psiCombustionModel/psiThermoCombustion/psiThermoCombustion.C
index abcb275c15..89a4a5daa7 100644
--- a/src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModel.C
+++ b/src/combustionModels/psiCombustionModel/psiThermoCombustion/psiThermoCombustion.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,54 +23,48 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "rhoChemistryCombustionModel.H"
-
-/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
-
-namespace Foam
-{
-namespace combustionModels
-{
-    defineTypeNameAndDebug(rhoChemistryCombustionModel, 0);
-    defineRunTimeSelectionTable(rhoChemistryCombustionModel, dictionary);
-
-
+#include "psiThermoCombustion.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-rhoChemistryCombustionModel::rhoChemistryCombustionModel
+Foam::combustionModels::psiThermoCombustion::psiThermoCombustion
 (
     const word& modelType,
     const fvMesh& mesh
 )
 :
-    combustionModel(modelType, mesh),
-    pChemistry_(rhoChemistryModel::New(mesh))
+    psiCombustionModel(modelType, mesh),
+    thermoPtr_(psiReactionThermo::New(mesh))
 {}
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-rhoChemistryCombustionModel::~rhoChemistryCombustionModel()
+Foam::combustionModels::psiThermoCombustion::~psiThermoCombustion()
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-
-bool rhoChemistryCombustionModel::read()
+Foam::psiReactionThermo&
+Foam::combustionModels::psiThermoCombustion::thermo()
 {
-    if (combustionModel::read())
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
+    return thermoPtr_();
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-} // End namespace combustionModels
+const Foam::psiReactionThermo&
+Foam::combustionModels::psiThermoCombustion::thermo() const
+{
+    return thermoPtr_();
+}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::psiThermoCombustion::rho() const
+{
+    return thermoPtr_->rho();
+}
+
+
+// ************************************************************************* //
diff --git a/src/combustionModels/psiCombustionModel/psiThermoCombustion/psiThermoCombustion.H b/src/combustionModels/psiCombustionModel/psiThermoCombustion/psiThermoCombustion.H
new file mode 100644
index 0000000000..8a294718a3
--- /dev/null
+++ b/src/combustionModels/psiCombustionModel/psiThermoCombustion/psiThermoCombustion.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::psiThermoCombustion
+
+Description
+    Compressibility-based thermo model wrapper for combustion models
+
+SourceFiles
+    psiThermoCombustion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef psiThermoCombustion_H
+#define psiThermoCombustion_H
+
+#include "autoPtr.H"
+#include "psiCombustionModel.H"
+#include "psiReactionThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    class psiThermoCombustion Declaration
+\*---------------------------------------------------------------------------*/
+
+class psiThermoCombustion
+:
+    public psiCombustionModel
+{
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        psiThermoCombustion(const psiThermoCombustion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const psiThermoCombustion&);
+
+
+protected:
+
+    // Protected data
+
+        //- Pointer to chemistry model
+        autoPtr<psiReactionThermo> thermoPtr_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components and thermo
+        psiThermoCombustion(const word& modelType, const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~psiThermoCombustion();
+
+
+    // Member Functions
+
+        //- Return access to the thermo package
+        virtual psiReactionThermo& thermo();
+
+        //- Return const access to the thermo package
+        virtual const psiReactionThermo& thermo() const;
+
+        //- Return const access to the density field
+        virtual tmp<volScalarField> rho() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModelNew.C b/src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModelNew.C
deleted file mode 100644
index 6bf82ede51..0000000000
--- a/src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModelNew.C
+++ /dev/null
@@ -1,78 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "rhoChemistryCombustionModel.H"
-
-// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::combustionModels::rhoChemistryCombustionModel>
-Foam::combustionModels::rhoChemistryCombustionModel::New
-(
-    const fvMesh& mesh
-)
-{
-    const word combModelName
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "combustionProperties",
-                mesh.time().constant(),
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        ).lookup("combustionModel")
-    );
-
-    Info<< "Selecting combustion model " << combModelName << endl;
-
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(combModelName);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
-    {
-        FatalErrorIn
-        (
-            "rhoChemistryCombustionModel::New"
-        )   << "Unknown rhoChemistryCombustionModel type "
-            << combModelName << endl << endl
-            << "Valid  combustionModels are : " << endl
-            << dictionaryConstructorTablePtr_->toc()
-            << exit(FatalError);
-    }
-
-    const label tempOpen = combModelName.find('<');
-
-    const word className = combModelName(0, tempOpen);
-
-    return autoPtr<rhoChemistryCombustionModel>
-        (cstrIter()(className, mesh));
-}
-
-
-// ************************************************************************* //
diff --git a/src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModelI.H b/src/combustionModels/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
similarity index 54%
rename from src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModelI.H
rename to src/combustionModels/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
index 4551d2e108..22495abb4a 100644
--- a/src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModelI.H
+++ b/src/combustionModels/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,37 +23,48 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+#include "rhoChemistryCombustion.H"
 
-inline Foam::rhoChemistryModel&
-Foam::combustionModels::rhoChemistryCombustionModel::pChemistry()
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoChemistryCombustion::rhoChemistryCombustion
+(
+    const word& modelType,
+    const fvMesh& mesh
+)
+:
+    rhoCombustionModel(modelType, mesh),
+    chemistryPtr_(rhoChemistryModel::New(mesh))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoChemistryCombustion::~rhoChemistryCombustion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::rhoReactionThermo&
+Foam::combustionModels::rhoChemistryCombustion::thermo()
 {
-    return pChemistry_();
+    return chemistryPtr_->thermo();
 }
 
-inline const Foam::rhoChemistryModel&
-Foam::combustionModels::rhoChemistryCombustionModel::
-pChemistry() const
+
+const Foam::rhoReactionThermo&
+Foam::combustionModels::rhoChemistryCombustion::thermo() const
 {
-    return pChemistry_();
+    return chemistryPtr_->thermo();
 }
 
-inline Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::rhoChemistryCombustionModel::rho() const
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::rhoChemistryCombustion::rho() const
 {
-    return pChemistry_->thermo().rho();
+    return chemistryPtr_->thermo().rho();
 }
 
-inline const Foam::rhoReactionThermo&
-Foam::combustionModels::rhoChemistryCombustionModel::thermo() const
-{
-    return pChemistry_->thermo();
-}
-
-inline Foam::rhoReactionThermo&
-Foam::combustionModels::rhoChemistryCombustionModel::thermo()
-{
-    return pChemistry_->thermo();
-}
 
 // ************************************************************************* //
diff --git a/src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModel.H b/src/combustionModels/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.H
similarity index 51%
rename from src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModel.H
rename to src/combustionModels/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.H
index 9a7aa0af3d..7010cb5389 100644
--- a/src/combustionModels/rhoChemistryCombustionModel/rhoChemistryCombustionModel.H
+++ b/src/combustionModels/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,24 +22,21 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::rhoChemistryCombustionModel
+    Foam::rhoChemistryCombustion
 
 Description
-    Combustion models for compressibility-based thermodynamics
+    Density-based chemistry model wrapper for combustion models
 
 SourceFiles
-    rhoChemistryCombustionModelI.H
-    rhoChemistryCombustionModel.C
-    rhoChemistryCombustionModelNew.C
+    rhoChemistryCombustion.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef rhoChemistryCombustionModel_H
-#define rhoChemistryCombustionModel_H
+#ifndef rhoChemistryCombustion_H
+#define rhoChemistryCombustion_H
 
-#include "combustionModel.H"
 #include "autoPtr.H"
-#include "runTimeSelectionTables.H"
+#include "rhoCombustionModel.H"
 #include "rhoChemistryModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -50,108 +47,59 @@ namespace combustionModels
 {
 
 /*---------------------------------------------------------------------------*\
-                     class rhoChemistryCombustionModel Declaration
+                    class rhoChemistryCombustion Declaration
 \*---------------------------------------------------------------------------*/
 
-class rhoChemistryCombustionModel
+class rhoChemistryCombustion
 :
-    public combustionModel
+    public rhoCombustionModel
 {
-
     // Private Member Functions
 
         //- Construct as copy (not implemented)
-        rhoChemistryCombustionModel(const rhoChemistryCombustionModel&);
+        rhoChemistryCombustion(const rhoChemistryCombustion&);
 
         //- Disallow default bitwise assignment
-        void operator=(const rhoChemistryCombustionModel&);
+        void operator=(const rhoChemistryCombustion&);
 
 
 protected:
 
     // Protected data
 
-        //- Auto pointer to psiChemistry
-        autoPtr<rhoChemistryModel> pChemistry_;
+        //- Pointer to chemistry model
+        autoPtr<rhoChemistryModel> chemistryPtr_;
 
 
 public:
 
-    //- Runtime type information
-    TypeName("rhoChemistryCombustionModel");
-
-
-    //- Declare run-time constructor selection tables
-    declareRunTimeSelectionTable
-    (
-        autoPtr,
-        rhoChemistryCombustionModel,
-        dictionary,
-        (
-            const word& modelType,
-            const fvMesh& mesh
-        ),
-        (modelType, mesh)
-    );
-
-
     // Constructors
 
-
         //- Construct from components and thermo
-        rhoChemistryCombustionModel
-        (
-            const word& modelType,
-            const fvMesh& mesh
-        );
-
-
-
-    //- Selector
-    static autoPtr<rhoChemistryCombustionModel> New
-    (
-        const fvMesh& mesh
-    );
+        rhoChemistryCombustion(const word& modelType, const fvMesh& mesh);
 
 
     //- Destructor
-    virtual ~rhoChemistryCombustionModel();
+    virtual ~rhoChemistryCombustion();
 
 
     // Member Functions
 
-
         //- Return access to the thermo package
-        inline rhoChemistryModel& pChemistry();
+        virtual rhoReactionThermo& thermo();
 
         //- Return const access to the thermo package
-        inline const rhoChemistryModel& pChemistry() const;
+        virtual const rhoReactionThermo& thermo() const;
 
-        //- Return const access to rho
-        inline tmp<volScalarField> rho() const;
-
-        //- Return const access to rho
-        inline const rhoReactionThermo& thermo() const;
-
-        //- Return non const access to rho
-        inline rhoReactionThermo& thermo();
-
-
-    // I-O
-
-        //- Update properties from given dictionary
-        virtual bool read();
+        //- Return const access to the density field
+        virtual tmp<volScalarField> rho() const;
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
 } // End namespace combustionModels
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "rhoChemistryCombustionModelI.H"
+} // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/combustionModels/rhoCombustionModel/rhoCombustionModel.C b/src/combustionModels/rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.C
similarity index 91%
rename from src/combustionModels/rhoCombustionModel/rhoCombustionModel.C
rename to src/combustionModels/rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.C
index 7a0a83a64d..55655caeeb 100644
--- a/src/combustionModels/rhoCombustionModel/rhoCombustionModel.C
+++ b/src/combustionModels/rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.C
@@ -45,24 +45,18 @@ Foam::combustionModels::rhoCombustionModel::rhoCombustionModel
     const fvMesh& mesh
 )
 :
-    combustionModel(modelType, mesh),
-    thermo_(rhoReactionThermo::New(mesh))
+    combustionModel(modelType, mesh)
 {}
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::combustionModels::rhoCombustionModel::~rhoCombustionModel()
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::rhoCombustionModel::rho() const
-{
-    return thermo_->rho();
-}
-
-
 bool Foam::combustionModels::rhoCombustionModel::read()
 {
     if (combustionModel::read())
@@ -75,4 +69,5 @@ bool Foam::combustionModels::rhoCombustionModel::read()
     }
 }
 
+
 // ************************************************************************* //
diff --git a/src/combustionModels/rhoCombustionModel/rhoCombustionModel.H b/src/combustionModels/rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.H
similarity index 85%
rename from src/combustionModels/rhoCombustionModel/rhoCombustionModel.H
rename to src/combustionModels/rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.H
index 8ab5e25740..04cf1f56ab 100644
--- a/src/combustionModels/rhoCombustionModel/rhoCombustionModel.H
+++ b/src/combustionModels/rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.H
@@ -67,14 +67,6 @@ class rhoCombustionModel
         void operator=(const rhoCombustionModel&);
 
 
-protected:
-
-    // Protected data
-
-        //- Thermo package
-        autoPtr<rhoReactionThermo> thermo_;
-
-
 public:
 
     //- Runtime type information
@@ -99,19 +91,12 @@ public:
 
 
         //- Construct from components and thermo
-        rhoCombustionModel
-        (
-            const word& modelType,
-            const fvMesh& mesh
-        );
+        rhoCombustionModel(const word& modelType, const fvMesh& mesh);
 
 
 
     //- Selector
-    static autoPtr<rhoCombustionModel> New
-    (
-        const fvMesh& mesh
-    );
+    static autoPtr<rhoCombustionModel> New(const fvMesh& mesh);
 
 
     //- Destructor
@@ -127,13 +112,14 @@ public:
 
 
         //- Return access to the thermo package
-        inline rhoReactionThermo& thermo();
+        virtual rhoReactionThermo& thermo() = 0;
 
         //- Return const access to the thermo package
-        inline const rhoReactionThermo& thermo() const;
+        virtual const rhoReactionThermo& thermo() const = 0;
 
         //- Return tmp of rho
-        virtual tmp<volScalarField> rho() const;
+        virtual tmp<volScalarField> rho() const = 0;
+
 
      // I-O
 
@@ -144,12 +130,8 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
 } // End namespace combustionModels
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "rhoCombustionModelI.H"
+} // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/combustionModels/rhoCombustionModel/rhoCombustionModelNew.C b/src/combustionModels/rhoCombustionModel/rhoCombustionModel/rhoCombustionModelNew.C
similarity index 100%
rename from src/combustionModels/rhoCombustionModel/rhoCombustionModelNew.C
rename to src/combustionModels/rhoCombustionModel/rhoCombustionModel/rhoCombustionModelNew.C
diff --git a/src/combustionModels/rhoCombustionModel/rhoCombustionModelI.H b/src/combustionModels/rhoCombustionModel/rhoCombustionModelI.H
deleted file mode 100644
index cbf0a65693..0000000000
--- a/src/combustionModels/rhoCombustionModel/rhoCombustionModelI.H
+++ /dev/null
@@ -1,41 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline Foam::rhoReactionThermo&
-Foam::combustionModels::rhoCombustionModel::thermo()
-{
-    return thermo_();
-}
-
-inline const Foam::rhoReactionThermo&
-Foam::combustionModels::rhoCombustionModel::thermo() const
-{
-    return thermo_();
-}
-
-
-// ************************************************************************* //
diff --git a/src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModelI.H b/src/combustionModels/rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.C
similarity index 55%
rename from src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModelI.H
rename to src/combustionModels/rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.C
index 30d1801739..da06dfda61 100644
--- a/src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModelI.H
+++ b/src/combustionModels/rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,37 +23,48 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+#include "rhoThermoCombustion.H"
 
-inline Foam::psiChemistryModel&
-Foam::combustionModels::psiChemistryCombustionModel::pChemistry()
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoThermoCombustion::rhoThermoCombustion
+(
+    const word& modelType,
+    const fvMesh& mesh
+)
+:
+    rhoCombustionModel(modelType, mesh),
+    thermoPtr_(rhoReactionThermo::New(mesh))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoThermoCombustion::~rhoThermoCombustion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::rhoReactionThermo&
+Foam::combustionModels::rhoThermoCombustion::thermo()
 {
-    return pChemistry_();
+    return thermoPtr_();
 }
 
-inline const Foam::psiChemistryModel&
-Foam::combustionModels::psiChemistryCombustionModel::
-pChemistry() const
+
+const Foam::rhoReactionThermo&
+Foam::combustionModels::rhoThermoCombustion::thermo() const
 {
-    return pChemistry_();
+    return thermoPtr_();
 }
 
-inline Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::psiChemistryCombustionModel::rho() const
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::rhoThermoCombustion::rho() const
 {
-    return pChemistry_->thermo().rho();
+    return thermoPtr_().rho();
 }
 
-inline const Foam::psiReactionThermo&
-Foam::combustionModels::psiChemistryCombustionModel::thermo() const
-{
-    return pChemistry_->thermo();
-}
-
-inline Foam::psiReactionThermo&
-Foam::combustionModels::psiChemistryCombustionModel::thermo()
-{
-    return pChemistry_->thermo();
-}
 
 // ************************************************************************* //
diff --git a/src/combustionModels/rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.H b/src/combustionModels/rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.H
new file mode 100644
index 0000000000..1a14440606
--- /dev/null
+++ b/src/combustionModels/rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::rhoThermoCombustion
+
+Description
+    Density-based thermo model wrapper for combustion models
+
+SourceFiles
+    rhoThermoCombustion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rhoThermoCombustion_H
+#define rhoThermoCombustion_H
+
+#include "autoPtr.H"
+#include "rhoCombustionModel.H"
+#include "rhoChemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    class rhoThermoCombustion Declaration
+\*---------------------------------------------------------------------------*/
+
+class rhoThermoCombustion
+:
+    public rhoCombustionModel
+{
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        rhoThermoCombustion(const rhoThermoCombustion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const rhoThermoCombustion&);
+
+
+protected:
+
+    // Protected data
+
+        //- Pointer to thermo model
+        autoPtr<rhoReactionThermo> thermoPtr_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components and thermo
+        rhoThermoCombustion(const word& modelType, const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~rhoThermoCombustion();
+
+
+    // Member Functions
+
+        //- Return access to the thermo package
+        virtual rhoReactionThermo& thermo();
+
+        //- Return const access to the thermo package
+        virtual const rhoReactionThermo& thermo() const;
+
+        //- Return const access to the density field
+        virtual tmp<volScalarField> rho() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/combustionModels/singleStepCombustion/singleStepCombustion.C b/src/combustionModels/singleStepCombustion/singleStepCombustion.C
index 319e5476b2..aab73f7abf 100644
--- a/src/combustionModels/singleStepCombustion/singleStepCombustion.C
+++ b/src/combustionModels/singleStepCombustion/singleStepCombustion.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,17 +33,14 @@ namespace combustionModels
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CombThermoType, class ThermoType>
-singleStepCombustion<CombThermoType, ThermoType>
-::singleStepCombustion
+singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
 (
-    const word& modelType, const fvMesh& mesh
+    const word& modelType,
+    const fvMesh& mesh
 )
 :
     CombThermoType(modelType, mesh),
-    singleMixture_
-    (
-        dynamic_cast<singleStepReactingMixture<ThermoType>&>(this->thermo())
-    ),
+    singleMixturePtr_(NULL),
     wFuel_
     (
         IOobject
@@ -55,22 +52,45 @@ singleStepCombustion<CombThermoType, ThermoType>
             IOobject::NO_WRITE
         ),
         this->mesh(),
-        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0)
+        dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
     )
-{}
+{
+    if (isA<singleStepReactingMixture<ThermoType> >(this->thermo()))
+    {
+        singleMixturePtr_ =
+            &dynamic_cast<singleStepReactingMixture<ThermoType>&>
+            (
+                this->thermo()
+            );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "singleStepCombustion<CombThermoType, ThermoType>::"
+            "singleStepCombustion"
+            "("
+                "const word&, "
+                "const fvMesh&"
+            ")"
+        )
+            << "Inconsistent thermo package for " << this->type() << " model:\n"
+            << "    " << this->thermo().type() << nl << nl
+            << "Please select a thermo package based on "
+            << "singleStepReactingMixture" << exit(FatalError);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
 
 template<class CombThermoType, class ThermoType>
-singleStepCombustion<CombThermoType, ThermoType>
-::~singleStepCombustion()
+singleStepCombustion<CombThermoType, ThermoType>::~singleStepCombustion()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-
 template<class CombThermoType, class ThermoType>
 Foam::tmp<Foam::fvScalarMatrix>
 singleStepCombustion<CombThermoType, ThermoType>::R
@@ -78,11 +98,11 @@ singleStepCombustion<CombThermoType, ThermoType>::R
     const volScalarField& Y
 ) const
 {
-    const label specieI = this->thermo_->composition().species()[Y.name()];
+    const label specieI = this->thermoPtr_->composition().species()[Y.name()];
 
     const volScalarField wSpecie
     (
-        wFuel_*singleMixture_.specieStoichCoeffs()[specieI]
+        wFuel_*singleMixturePtr_->specieStoichCoeffs()[specieI]
     );
 
     return wSpecie + fvm::Sp(0.0*wSpecie, Y);
@@ -93,10 +113,10 @@ template<class CombThermoType, class ThermoType>
 Foam::tmp<Foam::volScalarField>
 singleStepCombustion< CombThermoType, ThermoType>::Sh() const
 {
-    const label fuelI = singleMixture_.fuelIndex();
-    const volScalarField& YFuel = this->thermo_->composition().Y(fuelI);
+    const label fuelI = singleMixturePtr_->fuelIndex();
+    const volScalarField& YFuel = this->thermoPtr_->composition().Y(fuelI);
 
-    return -singleMixture_.qFuel()*(R(YFuel) & YFuel);
+    return -singleMixturePtr_->qFuel()*(R(YFuel) & YFuel);
 }
 
 
diff --git a/src/combustionModels/singleStepCombustion/singleStepCombustion.H b/src/combustionModels/singleStepCombustion/singleStepCombustion.H
index 0d50dfaefe..fd743247b3 100644
--- a/src/combustionModels/singleStepCombustion/singleStepCombustion.H
+++ b/src/combustionModels/singleStepCombustion/singleStepCombustion.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,7 +44,7 @@ namespace combustionModels
 {
 
 /*---------------------------------------------------------------------------*\
-                  Class singleStepCombustion Declaration
+                    Class singleStepCombustion Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CombThermoType, class ThermoType>
@@ -52,20 +52,6 @@ class singleStepCombustion
 :
     public CombThermoType
 {
-
-protected:
-
-    // Protected data
-
-        //- Reference to singleStepReactingMixture mixture
-        singleStepReactingMixture<ThermoType>& singleMixture_;
-
-        //- Fuel consumption rate
-        volScalarField wFuel_;
-
-
-private:
-
     // Private Member Functions
 
         //- Disallow copy construct
@@ -75,16 +61,23 @@ private:
         void operator=(const singleStepCombustion&);
 
 
-public:
+protected:
 
+    // Protected data
+
+        //- Pointer to singleStepReactingMixture mixture
+        singleStepReactingMixture<ThermoType>* singleMixturePtr_;
+
+        //- Fuel consumption rate
+        volScalarField wFuel_;
+
+
+public:
 
     // Constructors
 
         //- Construct from components
-        singleStepCombustion
-        (
-            const word& modelType, const fvMesh& mesh
-        );
+        singleStepCombustion(const word& modelType, const fvMesh& mesh);
 
 
     //- Destructor
@@ -95,7 +88,6 @@ public:
 
         // Evolution
 
-
             //- Fuel consumption rate matrix
             virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
 
@@ -106,7 +98,7 @@ public:
             virtual tmp<volScalarField> Sh() const;
 
 
-    // I-O
+        // I-O
 
             //- Update properties from given dictionary
             virtual bool read();
@@ -120,7 +112,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-
 #ifdef NoRepository
 #   include "singleStepCombustion.C"
 #endif
diff --git a/src/dynamicMesh/boundaryMesh/boundaryMesh.C b/src/dynamicMesh/boundaryMesh/boundaryMesh.C
index 51be127dd7..4a705acead 100644
--- a/src/dynamicMesh/boundaryMesh/boundaryMesh.C
+++ b/src/dynamicMesh/boundaryMesh/boundaryMesh.C
@@ -924,10 +924,10 @@ Foam::labelList Foam::boundaryMesh::getNearest
     // Create the octrees
     indexedOctree
     <
-        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+        treeDataPrimitivePatch<uindirectPrimitivePatch>
     > leftTree
     (
-        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+        treeDataPrimitivePatch<uindirectPrimitivePatch>
         (
             false,          // cacheBb
             leftPatch
@@ -939,10 +939,10 @@ Foam::labelList Foam::boundaryMesh::getNearest
     );
     indexedOctree
     <
-        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+        treeDataPrimitivePatch<uindirectPrimitivePatch>
     > rightTree
     (
-        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+        treeDataPrimitivePatch<uindirectPrimitivePatch>
         (
             false,          // cacheBb
             rightPatch
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
index c9d50b8930..56a81e0f61 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
@@ -209,8 +209,8 @@ void Foam::waveSurfacePressureFvPatchScalarField::updateCoeffs()
             << exit(FatalError);
     }
 
-    Info<< "min/max mag(zetap) = " << min(zetap & nf()) << ", "
-        << max(zetap & nf()) << endl;
+    Info<< "min/max(zetap) = " << gMin(zetap & nf()) << ", "
+        << gMax(zetap & nf()) << endl;
 
     fixedValueFvPatchScalarField::updateCoeffs();
 }
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateSurfaceInterpolationScheme/multivariateSurfaceInterpolationScheme.H b/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateSurfaceInterpolationScheme/multivariateSurfaceInterpolationScheme.H
index 00bc2b8877..aaea76b2a9 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateSurfaceInterpolationScheme/multivariateSurfaceInterpolationScheme.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateSurfaceInterpolationScheme/multivariateSurfaceInterpolationScheme.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,7 +44,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                 Class multivariateSurfaceInterpolationScheme Declaration
+           Class multivariateSurfaceInterpolationScheme Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index ea4b85b865..b64b1e8f99 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -318,7 +318,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     Sph = dt*htc*As;
 
-    dhsTrans += Sph*(0.5*(T_ + Tnew) - Tc_);
+    dhsTrans += Sph*(Tres.average() - Tc_);
 
     return Tnew;
 }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index a532285716..1c10870726 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -26,6 +26,7 @@ License
 #include "InjectionModel.H"
 #include "mathematicalConstants.H"
 #include "meshTools.H"
+#include "volFields.H"
 
 using namespace Foam::constant::mathematical;
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
index d626a9ae25..30a2c97b16 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
@@ -52,6 +52,7 @@ SourceFiles
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
 #include "SubModelBase.H"
+#include "vector.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
index faafdd551d..76a22ffe44 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
@@ -396,7 +396,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::H
             forAll(Y, i)
             {
                 label gid = props.globalIds()[i];
-                HMixture += Y[i]*thermo_.carrier().Hs(gid, p, T);
+                HMixture += Y[i]*thermo_.carrier().Ha(gid, p, T);
             }
             break;
         }
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 2e908b7318..35d0822d2a 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -249,9 +249,10 @@ class autoSnapDriver
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
 
-                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceSurfNormals,
                     const List<List<point> >& pointFaceDisp,
                     const List<List<point> >& pointFaceCentres,
+                    const labelListList& pointFacePatchID,
 
                     vectorField& patchAttraction,
                     List<pointConstraint>& patchConstraints
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 6347d5dd16..950849e15e 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -71,14 +71,15 @@ namespace Foam
         }
     };
 
+    template<class T>
     class listPlusEqOp
     {
     public:
 
         void operator()
         (
-            List<point>& x,
-            const List<point>& y
+            List<T>& x,
+            const List<T>& y
         ) const
         {
             label sz = x.size();
@@ -486,7 +487,7 @@ void Foam::autoSnapDriver::binFeatureFaces
 
     const label pointI,
 
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
 
@@ -495,7 +496,7 @@ void Foam::autoSnapDriver::binFeatureFaces
     DynamicList<label>& surfaceCount
 ) const
 {
-    const List<point>& pfNormals = pointFaceNormals[pointI];
+    const List<point>& pfNormals = pointFaceSurfNormals[pointI];
     const List<point>& pfDisp = pointFaceDisp[pointI];
     const List<point>& pfCentres = pointFaceCentres[pointI];
 
@@ -531,7 +532,7 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
 
@@ -590,7 +591,7 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
             snapDist,
             pointI,
 
-            pointFaceNormals,
+            pointFaceSurfNormals,
             pointFaceDisp,
             pointFaceCentres,
 
@@ -720,7 +721,7 @@ void Foam::autoSnapDriver::determineAllFeatures
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
 
@@ -873,7 +874,7 @@ void Foam::autoSnapDriver::determineFeatures
     //const vectorField& faceSurfaceNormal,
     //const vectorField& faceDisp,
     //const vectorField& faceRotation,
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
 
@@ -955,7 +956,7 @@ void Foam::autoSnapDriver::determineFeatures
 
                 //faceSurfaceNormal,
                 //faceDisp,
-                pointFaceNormals,
+                pointFaceSurfNormals,
                 pointFaceDisp,
                 pointFaceCentres,
 
@@ -1132,9 +1133,10 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     //const vectorField& faceSurfaceNormal,
     //const vectorField& faceDisp,
     //const vectorField& faceRotation,
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
+    const labelListList& pointFacePatchID,
 
     vectorField& patchAttraction,
     List<pointConstraint>& patchConstraints
@@ -1181,7 +1183,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
         pp,
         snapDist,
 
-        pointFaceNormals,
+        pointFaceSurfNormals,
         pointFaceDisp,
         pointFaceCentres,
 
@@ -1340,6 +1342,131 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     }
 
 
+
+    //MEJ: any faces that have multi-patch points only keep the multi-patch
+    //     points.
+    {
+        autoPtr<OFstream> multiPatchStr;
+        if (debug&meshRefinement::OBJINTERSECTIONS)
+        {
+            multiPatchStr.reset
+            (
+                new OFstream
+                (
+                    meshRefiner_.mesh().time().path()
+                  / "multiPatch_" + name(iter) + ".obj"
+                )
+            );
+            Pout<< nl << "Dumping removed constraints due to same-face"
+                << " multi-patch points to "
+                << multiPatchStr().name() << endl;
+        }
+
+
+        // 1. Mark points on multiple patches
+        PackedBoolList isMultiPatchPoint(pp.size());
+
+        forAll(pointFacePatchID, pointI)
+        {
+            const labelList& patches = pointFacePatchID[pointI];
+            label patch0 = patches[0];
+            for (label i = 1; i < patches.size(); i++)
+            {
+                if (patch0 != patches[i])
+                {
+                    isMultiPatchPoint[pointI] = 1u;
+                    break;
+                }
+            }
+        }
+
+        // 2. Make sure multi-patch points are also attracted
+        forAll(isMultiPatchPoint, pointI)
+        {
+            if (isMultiPatchPoint[pointI])
+            {
+                if
+                (
+                    patchConstraints[pointI].first() == 0
+                 && allPatchConstraints[pointI].first() > 0
+                )
+                {
+                    patchAttraction[pointI] = allPatchAttraction[pointI];
+                    patchConstraints[pointI] = allPatchConstraints[pointI];
+
+                    if (multiPatchStr.valid())
+                    {
+                        Pout<< "Adding constraint on multiPatchPoint:"
+                            << pp.localPoints()[pointI]
+                            << " constraint:" << patchConstraints[pointI]
+                            << " attraction:" << patchAttraction[pointI]
+                            << endl;
+                    }
+                }
+            }
+        }
+
+        // Up to here it is all parallel ok.
+
+
+        // 3. Knock out any attraction on faces with multi-patch points
+        label nChanged = 0;
+        forAll(pp.localFaces(), faceI)
+        {
+            const face& f = pp.localFaces()[faceI];
+
+            label nMultiPatchPoints = 0;
+            forAll(f, fp)
+            {
+                label pointI = f[fp];
+                if
+                (
+                    isMultiPatchPoint[pointI]
+                 && patchConstraints[pointI].first() != 0
+                )
+                {
+                    nMultiPatchPoints++;
+                }
+            }
+
+            if (nMultiPatchPoints > 0)
+            {
+                forAll(f, fp)
+                {
+                    label pointI = f[fp];
+                    if
+                    (
+                       !isMultiPatchPoint[pointI]
+                     && patchConstraints[pointI].first() != 0
+                    )
+                    {
+                        //Pout<< "Knocking out constraint"
+                        //    << " on non-multiPatchPoint:"
+                        //    << pp.localPoints()[pointI] << endl;
+                        patchAttraction[pointI] = vector::zero;
+                        patchConstraints[pointI] = pointConstraint();
+                        nChanged++;
+
+                        if (multiPatchStr.valid())
+                        {
+                            meshTools::writeOBJ
+                            (
+                                multiPatchStr(),
+                                pp.localPoints()[pointI]
+                            );
+                        }
+                    }
+                }
+            }
+        }
+
+        reduce(nChanged, sumOp<label>());
+        Info<< "Removing constraints near multi-patch points : changed "
+            << nChanged << " points" << endl;
+    }
+
+
+
     // Dump
     if (debug&meshRefinement::OBJINTERSECTIONS)
     {
@@ -1753,8 +1880,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     const pointField& localPoints = pp.localPoints();
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    // Displacement and orientation per pp face.
+    // Displacement and orientation per pp face
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // vector from point on surface back to face centre
     vectorField faceDisp(pp.size(), vector::zero);
+    // normal of surface at point on surface
     vectorField faceSurfaceNormal(pp.size(), vector::zero);
     vectorField faceRotation(pp.size(), vector::zero);
 
@@ -1771,34 +1902,42 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     vectorField patchDisp = nearestDisp;
 
 
-    // Collect (possibly remote) face-wise data on coupled points.
+    // Collect (possibly remote) per point data of all surrounding faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // - faceSurfaceNormal
     // - faceDisp
     // - faceRotation
-    // - faceCentres
+    // - faceCentres&faceNormal
 
     // For now just get all surrounding face data. Expensive - should just
     // store and sync data on coupled points only (see e.g PatchToolsNormals.C)
 
-    List<List<point> > pointFaceNormals(pp.nPoints());
+    List<List<point> > pointFaceSurfNormals(pp.nPoints());
     List<List<point> > pointFaceDisp(pp.nPoints());
     List<List<point> > pointFaceCentres(pp.nPoints());
+    List<labelList>    pointFacePatchID(pp.nPoints());
 
     // Fill local data
     forAll(pp.pointFaces(), pointI)
     {
         const labelList& pFaces = pp.pointFaces()[pointI];
-        List<point>& pNormals = pointFaceNormals[pointI];
+        List<point>& pNormals = pointFaceSurfNormals[pointI];
         pNormals.setSize(pFaces.size());
         List<point>& pDisp = pointFaceDisp[pointI];
         pDisp.setSize(pFaces.size());
         List<point>& pFc = pointFaceCentres[pointI];
         pFc.setSize(pFaces.size());
+        labelList& pFid = pointFacePatchID[pointI];
+        pFid.setSize(pFaces.size());
+
         forAll(pFaces, i)
         {
-            pNormals[i] = faceSurfaceNormal[pFaces[i]];
-            pDisp[i] = faceDisp[pFaces[i]];
-            pFc[i] = pp.faceCentres()[pFaces[i]];
+            label faceI = pFaces[i];
+            pNormals[i] = faceSurfaceNormal[faceI];
+            pDisp[i] = faceDisp[faceI];
+            pFc[i] = pp.faceCentres()[faceI];
+            label meshFaceI = pp.addressing()[faceI];
+            pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
         }
     }
 
@@ -1806,8 +1945,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     (
         mesh,
         pp.meshPoints(),
-        pointFaceNormals,
-        listPlusEqOp(),
+        pointFaceSurfNormals,
+        listPlusEqOp<point>(),
         List<point>(),
         listTransform()
     );
@@ -1816,7 +1955,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         mesh,
         pp.meshPoints(),
         pointFaceDisp,
-        listPlusEqOp(),
+        listPlusEqOp<point>(),
         List<point>(),
         listTransform()
     );
@@ -1825,13 +1964,27 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         mesh,
         pp.meshPoints(),
         pointFaceCentres,
-        listPlusEqOp(),
+        listPlusEqOp<point>(),
         List<point>(),
         listTransform()
     );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFacePatchID,
+        listPlusEqOp<label>(),
+        List<label>()
+    );
 
 
 
+    // Main calculation
+    // ~~~~~~~~~~~~~~~~
+    // This is the main intelligence which calculates per point the vector to
+    // attract it to the nearest surface. There are lots of possibilities
+    // here.
+
     // Nearest feature
     vectorField patchAttraction(localPoints.size(), vector::zero);
     // Constraints at feature
@@ -1845,9 +1998,10 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         pp,
         snapDist,
 
-        pointFaceNormals,
+        pointFaceSurfNormals,
         pointFaceDisp,
         pointFaceCentres,
+        pointFacePatchID,
 
         patchAttraction,
         patchConstraints
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 2664719b26..e1b76c44e1 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -272,7 +272,11 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
     // Find all start cells of features. Is done by tracking from keepPoint.
     Cloud<trackedParticle> cloud(mesh_, IDLList<trackedParticle>());
 
-    // Create particles on whichever processor holds the keepPoint.
+
+    // Features are identical on all processors. Number them so we know
+    // what to seed. Do this on only the processor that
+    // holds the keepPoint.
+
     label cellI = -1;
     label tetFaceI = -1;
     label tetPtI = -1;
@@ -281,13 +285,25 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
 
     if (cellI != -1)
     {
+        // I am the processor that holds the keepPoint
+
         forAll(features_, featI)
         {
             const featureEdgeMesh& featureMesh = features_[featI];
             const label featureLevel = features_.levels()[featI];
-
             const labelListList& pointEdges = featureMesh.pointEdges();
 
+            // Find regions on edgeMesh
+            labelList edgeRegion;
+            label nRegions = featureMesh.regions(edgeRegion);
+
+
+            PackedBoolList regionVisited(nRegions);
+
+
+            // 1. Seed all 'knots' in edgeMesh
+
+
             forAll(pointEdges, pointI)
             {
                 if (pointEdges[pointI].size() != 2)
@@ -316,6 +332,50 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
                             pointI                          // end point
                         )
                     );
+
+                    // Mark
+                    if (pointEdges[pointI].size() > 0)
+                    {
+                        label e0 = pointEdges[pointI][0];
+                        label regionI = edgeRegion[e0];
+                        regionVisited[regionI] = 1u;
+                    }
+                }
+            }
+
+
+            // 2. Any regions that have not been visited at all? These can
+            //    only be circular regions!
+            forAll(featureMesh.edges(), edgeI)
+            {
+                if (regionVisited.set(edgeRegion[edgeI], 1u))
+                {
+                    const edge& e = featureMesh.edges()[edgeI];
+                    label pointI = e.start();
+                    if (debug)
+                    {
+                        Pout<< "Adding particle from point:" << pointI
+                            << " coord:" << featureMesh.points()[pointI]
+                            << " on circular region:" << edgeRegion[edgeI]
+                            << endl;
+                    }
+
+                    // Non-manifold point. Create particle.
+                    cloud.addParticle
+                    (
+                        new trackedParticle
+                        (
+                            mesh_,
+                            keepPoint,
+                            cellI,
+                            tetFaceI,
+                            tetPtI,
+                            featureMesh.points()[pointI],   // endpos
+                            featureLevel,                   // level
+                            featI,                          // featureMesh
+                            pointI                          // end point
+                        )
+                    );
                 }
             }
         }
@@ -329,6 +389,8 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
     trackedParticle::trackingData td(cloud, maxFeatureLevel);
 
     // Track all particles to their end position (= starting feature point)
+    // Note that the particle might have started on a different processor
+    // so this will transfer across nicely until we can start tracking proper.
     cloud.move(td, GREAT);
 
     // Reset level
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index 3962697845..ad5a974ed7 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -149,13 +149,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::checkPatches
     const scalar maxBoundsError = 0.05;
 
     // check bounds of source and target
-    boundBox bbSrc(srcPatch.points(), srcPatch.meshPoints());
-    reduce(bbSrc.min(), minOp<point>());
-    reduce(bbSrc.max(), maxOp<point>());
-
-    boundBox bbTgt(tgtPatch.points(), tgtPatch.meshPoints());
-    reduce(bbTgt.min(), minOp<point>());
-    reduce(bbTgt.max(), maxOp<point>());
+    boundBox bbSrc(srcPatch.points(), srcPatch.meshPoints(), true);
+    boundBox bbTgt(tgtPatch.points(), tgtPatch.meshPoints(), true);
 
     boundBox bbTgtInf(bbTgt);
     bbTgtInf.inflate(maxBoundsError);
@@ -1164,50 +1159,36 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
     const word& patchName,
     const labelListList& addr,
     scalarListList& wght,
+    scalarField& wghtSum,
     const bool output
 )
 {
-    scalar minBound = VGREAT;
-    scalar maxBound = -VGREAT;
-
-    scalar tSum = 0.0;
-
     // Normalise the weights
+    wghtSum.setSize(wght.size());
     forAll(wght, faceI)
     {
         scalar s = sum(wght[faceI]);
         scalar t = s/patchAreas[faceI];
 
-        tSum += t;
-
-        if (t < minBound)
-        {
-            minBound = t;
-        }
-
-        if (t > maxBound)
-        {
-            maxBound = t;
-        }
-
         forAll(addr[faceI], i)
         {
             wght[faceI][i] /= s;
         }
+
+        wghtSum[faceI] = t;
     }
 
 
     if (output)
     {
-        const label nFace = returnReduce(wght.size(), sumOp<scalar>());
-        reduce(tSum, sumOp<scalar>());
+        const label nFace = returnReduce(wght.size(), sumOp<label>());
 
         if (nFace)
         {
             Info<< "AMI: Patch " << patchName << " weights min/max/average = "
-                << returnReduce(minBound, minOp<scalar>()) << ", "
-                << returnReduce(maxBound, maxOp<scalar>()) << ", "
-                << tSum/nFace << endl;
+                << gMin(wghtSum) << ", "
+                << gMax(wghtSum) << ", "
+                << gAverage(wghtSum) << endl;
         }
     }
 }
@@ -1227,6 +1208,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
     scalarField& srcMagSf,
     labelListList& srcAddress,
     scalarListList& srcWeights,
+    scalarField& srcWeightsSum,
     autoPtr<mapDistribute>& tgtMap
 )
 {
@@ -1468,6 +1450,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
         "source",
         srcAddress,
         srcWeights,
+        srcWeightsSum,
         false
     );
 }
@@ -1488,9 +1471,11 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     singlePatchProc_(-999),
     srcAddress_(),
     srcWeights_(),
+    srcWeightsSum_(),
     srcNonOverlap_(),
     tgtAddress_(),
     tgtWeights_(),
+    tgtWeightsSum_(),
     treePtr_(NULL),
     triMode_(triMode),
     srcMapPtr_(NULL),
@@ -1521,9 +1506,11 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     singlePatchProc_(-999),
     srcAddress_(),
     srcWeights_(),
+    srcWeightsSum_(),
     srcNonOverlap_(),
     tgtAddress_(),
     tgtWeights_(),
+    tgtWeightsSum_(),
     treePtr_(NULL),
     triMode_(triMode),
     srcMapPtr_(NULL),
@@ -1609,9 +1596,11 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     singlePatchProc_(fineAMI.singlePatchProc_),
     srcAddress_(),
     srcWeights_(),
+    srcWeightsSum_(),
     srcNonOverlap_(),
     tgtAddress_(),
     tgtWeights_(),
+    tgtWeightsSum_(),
     treePtr_(NULL),
     triMode_(fineAMI.triMode_),
     srcMapPtr_(NULL),
@@ -1681,6 +1670,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
         srcMagSf_,
         srcAddress_,
         srcWeights_,
+        srcWeightsSum_,
         tgtMapPtr_
     );
 
@@ -1706,6 +1696,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
         tgtMagSf_,
         tgtAddress_,
         tgtWeights_,
+        tgtWeightsSum_,
         srcMapPtr_
     );
 
@@ -1857,8 +1848,24 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         );
 
         // weights normalisation
-        normaliseWeights(srcMagSf_, "source", srcAddress_, srcWeights_, true);
-        normaliseWeights(tgtMagSf_, "target", tgtAddress_, tgtWeights_, true);
+        normaliseWeights
+        (
+            srcMagSf_,
+            "source",
+            srcAddress_,
+            srcWeights_,
+            srcWeightsSum_,
+            true
+        );
+        normaliseWeights
+        (
+            tgtMagSf_,
+            "target",
+            tgtAddress_,
+            tgtWeights_,
+            tgtWeightsSum_,
+            true
+        );
 
         // cache maps and reset addresses
         List<Map<label> > cMap;
@@ -1877,8 +1884,24 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
 
         calcAddressing(srcPatch, tgtPatch);
 
-        normaliseWeights(srcMagSf_, "source", srcAddress_, srcWeights_, true);
-        normaliseWeights(tgtMagSf_, "target", tgtAddress_, tgtWeights_, true);
+        normaliseWeights
+        (
+            srcMagSf_,
+            "source",
+            srcAddress_,
+            srcWeights_,
+            srcWeightsSum_,
+            true
+        );
+        normaliseWeights
+        (
+            tgtMagSf_,
+            "target",
+            tgtAddress_,
+            tgtWeights_,
+            tgtWeightsSum_,
+            true
+        );
     }
 
     if (debug)
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
index 36cd661175..b393030823 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
@@ -82,7 +82,7 @@ class AMIInterpolation
     public AMIInterpolationName
 {
     //- local typedef to octree tree-type
-    typedef treeDataPrimitivePatch<face, SubList, const pointField&> treeType;
+    typedef treeDataPrimitivePatch<TargetPatch> treeType;
 
 
     // Private data
@@ -106,6 +106,9 @@ class AMIInterpolation
             //- Weights of target faces per source face
             scalarListList srcWeights_;
 
+            //- Sum of weights of target faces per source face
+            scalarField srcWeightsSum_;
+
             //- Labels of faces that are not overlapped by any target faces
             //  (should be empty for correct functioning)
             labelList srcNonOverlap_;
@@ -122,6 +125,9 @@ class AMIInterpolation
             //- Weights of source faces per target face
             scalarListList tgtWeights_;
 
+            //- Sum of weights of source faces per target face
+            scalarField tgtWeightsSum_;
+
 
         //- Octree used to find face seeds
         autoPtr<indexedOctree<treeType> > treePtr_;
@@ -285,6 +291,7 @@ class AMIInterpolation
                 const word& patchName,
                 const labelListList& addr,
                 scalarListList& wght,
+                scalarField& wghtSum,
                 const bool output
             );
 
@@ -304,6 +311,7 @@ class AMIInterpolation
                 scalarField& srcMagSf,
                 labelListList& srcAddress,
                 scalarListList& srcWeights,
+                scalarField& srcWeightsSum,
                 autoPtr<mapDistribute>& tgtMap
             );
 
@@ -365,11 +373,14 @@ public:
                 //- Return const access to source patch weights
                 inline const scalarListList& srcWeights() const;
 
+                //- Return const access to normalisation factor of source
+                //  patch weights (i.e. the sum before normalisation)
+                inline const scalarField& srcWeightsSum() const;
+
                 //- Labels of faces that are not overlapped by any target faces
                 //  (should be empty for correct functioning)
                 inline const labelList& srcNonOverlap() const;
 
-
                 //- Source map pointer - valid only if singlePatchProc = -1
                 //  This gets source data into a form to be consumed by
                 //  tgtAddress, tgtWeights
@@ -387,7 +398,11 @@ public:
                 //- Return const access to target patch weights
                 inline const scalarListList& tgtWeights() const;
 
-                //- Target map pointer - valid only if singlePatchProc = -1
+                //- Return const access to normalisation factor of target
+                //  patch weights (i.e. the sum before normalisation)
+                inline const scalarField& tgtWeightsSum() const;
+
+                //- Target map pointer -  valid only if singlePatchProc=-1.
                 //  This gets target data into a form to be consumed by
                 //  srcAddress, srcWeights
                 inline const mapDistribute& tgtMap() const;
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
index 3dd21b5aa3..f2782bad88 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
@@ -63,6 +63,14 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcNonOverlap() const
 }
 
 
+template<class SourcePatch, class TargetPatch>
+inline const Foam::scalarField&
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeightsSum() const
+{
+    return srcWeightsSum_;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 inline const Foam::mapDistribute&
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMap() const
@@ -95,6 +103,14 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtWeights() const
 }
 
 
+template<class SourcePatch, class TargetPatch>
+inline const Foam::scalarField&
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtWeightsSum() const
+{
+    return tgtWeightsSum_;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 inline const Foam::mapDistribute&
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtMap() const
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 3329c2b28d..4abc57669e 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -110,6 +110,8 @@ $(faceSources)/patchToFace/patchToFace.C
 $(faceSources)/boundaryToFace/boundaryToFace.C
 $(faceSources)/zoneToFace/zoneToFace.C
 $(faceSources)/boxToFace/boxToFace.C
+$(faceSources)/regionToFace/regionToFace.C
+$(faceSources)/regionToFace/patchEdgeFaceRegion.C
 
 pointSources = sets/pointSources
 $(pointSources)/labelToPoint/labelToPoint.C
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
index 6927cbaa50..4cd2968232 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
@@ -29,30 +29,14 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::scalar
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-tolSqr = sqr(1e-6);
+template<class PatchType>
+Foam::scalar Foam::treeDataPrimitivePatch<PatchType>::tolSqr = sqr(1e-6);
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::treeBoundBox
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-calcBb
+template<class PatchType>
+Foam::treeBoundBox Foam::treeDataPrimitivePatch<PatchType>::calcBb
 (
     const pointField& points,
     const face& f
@@ -71,15 +55,8 @@ calcBb
 }
 
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-void Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-update()
+template<class PatchType>
+void Foam::treeDataPrimitivePatch<PatchType>::update()
 {
     if (cacheBb_)
     {
@@ -96,18 +73,11 @@ update()
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-treeDataPrimitivePatch
+template<class PatchType>
+Foam::treeDataPrimitivePatch<PatchType>::treeDataPrimitivePatch
 (
     const bool cacheBb,
-    const PrimitivePatch<Face, FaceList, PointField, PointType>& patch
+    const PatchType& patch
 )
 :
     patch_(patch),
@@ -119,16 +89,8 @@ treeDataPrimitivePatch
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::pointField
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-shapePoints() const
+template<class PatchType>
+Foam::pointField Foam::treeDataPrimitivePatch<PatchType>::shapePoints() const
 {
     pointField cc(patch_.size());
 
@@ -143,27 +105,10 @@ shapePoints() const
 
 //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
 //  Only makes sense for closed surfaces.
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::label
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-getVolumeType
+template<class PatchType>
+Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
 (
-    const indexedOctree
-    <
-        treeDataPrimitivePatch
-        <
-            Face,
-            FaceList,
-            PointField,
-            PointType
-        >
-    >& oc,
+    const indexedOctree<treeDataPrimitivePatch<PatchType> >& oc,
     const point& sample
 ) const
 {
@@ -396,16 +341,8 @@ getVolumeType
 
 
 // Check if any point on shape is inside cubeBb.
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-bool
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-overlaps
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
 (
     const label index,
     const treeBoundBox& cubeBb
@@ -462,16 +399,8 @@ overlaps
 
 
 // Check if any point on shape is inside sphere.
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-bool
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-overlaps
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
 (
     const label index,
     const point& centre,
@@ -512,16 +441,8 @@ overlaps
 
 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
 // nearestPoint.
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-void
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-findNearest
+template<class PatchType>
+void Foam::treeDataPrimitivePatch<PatchType>::findNearest
 (
     const labelUList& indices,
     const point& sample,
@@ -552,16 +473,8 @@ findNearest
 }
 
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-bool
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-intersects
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::intersects
 (
     const label index,
     const point& start,
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
index 44f76d2c75..bbb3c74956 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,13 +59,7 @@ TemplateName(treeDataPrimitivePatch);
                    Class treeDataPrimitivePatch Declaration
 \*---------------------------------------------------------------------------*/
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType=point
->
+template<class PatchType>
 class treeDataPrimitivePatch
 :
     public treeDataPrimitivePatchName
@@ -78,7 +72,7 @@ class treeDataPrimitivePatch
     // Private data
 
         //- Underlying geometry
-        const PrimitivePatch<Face, FaceList, PointField, PointType>& patch_;
+        const PatchType& patch_;
 
         //- Whether to precalculate and store face bounding box
         const bool cacheBb_;
@@ -103,7 +97,7 @@ public:
         treeDataPrimitivePatch
         (
             const bool cacheBb,
-            const PrimitivePatch<Face, FaceList, PointField, PointType>&
+            const PatchType&
         );
 
 
@@ -121,8 +115,7 @@ public:
             pointField shapePoints() const;
 
             //- Return access to the underlying patch
-            const PrimitivePatch<Face, FaceList, PointField, PointType>&
-            patch() const
+            const PatchType& patch() const
             {
                 return patch_;
             }
@@ -134,16 +127,7 @@ public:
             //  Only makes sense for closed surfaces.
             label getVolumeType
             (
-                const indexedOctree
-                <
-                    treeDataPrimitivePatch
-                    <
-                        Face,
-                        FaceList,
-                        PointField,
-                        PointType
-                    >
-                >&,
+                const indexedOctree<treeDataPrimitivePatch<PatchType> >&,
                 const point&
             ) const;
 
diff --git a/src/meshTools/searchableSurface/searchableSurface.C b/src/meshTools/searchableSurface/searchableSurface.C
index b586ca191d..decd4c7259 100644
--- a/src/meshTools/searchableSurface/searchableSurface.C
+++ b/src/meshTools/searchableSurface/searchableSurface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,4 +95,21 @@ Foam::searchableSurface::~searchableSurface()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::searchableSurface::findNearest
+(
+    const pointField& sample,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info,
+    vectorField& normal,
+    labelList& region
+) const
+{
+    findNearest(sample, nearestDistSqr, info);
+    getNormal(info, normal);
+    getRegion(info, region);
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableSurface.H b/src/meshTools/searchableSurface/searchableSurface.H
index 535997e900..94713fbf38 100644
--- a/src/meshTools/searchableSurface/searchableSurface.H
+++ b/src/meshTools/searchableSurface/searchableSurface.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -323,6 +323,17 @@ public:
                 List<volumeType>&
             ) const = 0;
 
+            //- Find nearest, normal and region. Can be overridden with
+            //  optimised implementation
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&,
+                vectorField& normal,
+                labelList& region
+            ) const;
+
 
         // Other
 
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index 6bbfd22960..86fcd5889b 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -827,6 +827,8 @@ void Foam::triSurfaceMesh::getNormal
     vectorField& normal
 ) const
 {
+    const triSurface& s = static_cast<const triSurface&>(*this);
+
     normal.setSize(info.size());
 
     if (minQuality_ >= 0)
@@ -834,7 +836,6 @@ void Foam::triSurfaceMesh::getNormal
         // Make sure we don't use triangles with low quality since
         // normal is not reliable.
 
-        const triSurface& s = static_cast<const triSurface&>(*this);
         const labelListList& faceFaces = s.faceFaces();
 
         forAll(info, i)
@@ -842,6 +843,7 @@ void Foam::triSurfaceMesh::getNormal
             if (info[i].hit())
             {
                 label faceI = info[i].index();
+                normal[i] = s[faceI].normal(points());
 
                 scalar qual = s[faceI].tri(points()).quality();
 
@@ -861,11 +863,8 @@ void Foam::triSurfaceMesh::getNormal
                         }
                     }
                 }
-                else
-                {
-                    normal[i] = s[faceI].normal(points());
-                }
-                normal[i] /= mag(normal[i]);
+
+                normal[i] /= mag(normal[i]) + VSMALL;
             }
             else
             {
@@ -885,7 +884,7 @@ void Foam::triSurfaceMesh::getNormal
                 //normal[i] = faceNormals()[faceI];
 
                 //- Uncached
-                normal[i] = triSurface::operator[](faceI).normal(points());
+                normal[i] = s[faceI].normal(points());
                 normal[i] /= mag(normal[i]) + VSMALL;
             }
             else
diff --git a/src/combustionModels/psiCombustionModel/psiCombustionModelI.H b/src/meshTools/sets/faceSources/regionToFace/patchEdgeFaceRegion.C
similarity index 73%
rename from src/combustionModels/psiCombustionModel/psiCombustionModelI.H
rename to src/meshTools/sets/faceSources/regionToFace/patchEdgeFaceRegion.C
index 5df2365b9d..91e9b4c331 100644
--- a/src/combustionModels/psiCombustionModel/psiCombustionModelI.H
+++ b/src/meshTools/sets/faceSources/regionToFace/patchEdgeFaceRegion.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,18 +23,27 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+#include "patchEdgeFaceRegion.H"
 
-inline Foam::psiReactionThermo&
-Foam::combustionModels::psiCombustionModel::thermo()
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::patchEdgeFaceRegion& wDist
+)
 {
-    return thermo_();
+    return os << wDist.region_;
 }
 
-inline const Foam::psiReactionThermo&
-Foam::combustionModels::psiCombustionModel::thermo() const
+
+Foam::Istream& Foam::operator>>
+(
+    Foam::Istream& is,
+    Foam::patchEdgeFaceRegion& wDist
+)
 {
-    return thermo_();
+    return is >> wDist.region_;
 }
 
 
diff --git a/src/meshTools/sets/faceSources/regionToFace/patchEdgeFaceRegion.H b/src/meshTools/sets/faceSources/regionToFace/patchEdgeFaceRegion.H
new file mode 100644
index 0000000000..6c0b65f2d7
--- /dev/null
+++ b/src/meshTools/sets/faceSources/regionToFace/patchEdgeFaceRegion.H
@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::patchEdgeFaceRegion
+
+Description
+    Transport of region for use in PatchEdgeFaceWave.
+
+    Set element to -2 to denote blocked.
+
+SourceFiles
+    patchEdgeFaceRegionI.H
+    patchEdgeFaceRegion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchEdgeFaceRegion_H
+#define patchEdgeFaceRegion_H
+
+#include "point.H"
+#include "label.H"
+#include "scalar.H"
+#include "tensor.H"
+#include "indirectPrimitivePatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class patchEdgeFaceRegion Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchEdgeFaceRegion
+{
+    // Private data
+
+        //- region
+        label region_;
+
+    // Private Member Functions
+
+        //- Combine current with w2. Update region_ if w2 has smaller
+        //  quantities and returns true.
+        template<class TrackingData>
+        inline bool update
+        (
+            const patchEdgeFaceRegion& w2,
+            const scalar tol,
+            TrackingData& td
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline patchEdgeFaceRegion();
+
+        //- Construct from origin, distance
+        inline patchEdgeFaceRegion(const label);
+
+
+    // Member Functions
+
+        // Access
+
+            inline label region() const;
+
+
+        // Needed by meshWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value.
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
+
+            //- Apply rotation matrix
+            template<class TrackingData>
+            inline void transform
+            (
+                const polyMesh& mesh,
+                const indirectPrimitivePatch& patch,
+                const tensor& rotTensor,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of face on edge
+            template<class TrackingData>
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const indirectPrimitivePatch& patch,
+                const label edgeI,
+                const label faceI,
+                const patchEdgeFaceRegion& faceInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- New information for edge (from e.g. coupled edge)
+            template<class TrackingData>
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const indirectPrimitivePatch& patch,
+                const patchEdgeFaceRegion& edgeInfo,
+                const bool sameOrientation,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of edge on face.
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh& mesh,
+                const indirectPrimitivePatch& patch,
+                const label faceI,
+                const label edgeI,
+                const patchEdgeFaceRegion& edgeInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Same (like operator==)
+            template<class TrackingData>
+            inline bool equal(const patchEdgeFaceRegion&, TrackingData&) const;
+
+
+    // Member Operators
+
+        // Needed for List IO
+        inline bool operator==(const patchEdgeFaceRegion&) const;
+        inline bool operator!=(const patchEdgeFaceRegion&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const patchEdgeFaceRegion&);
+        friend Istream& operator>>(Istream&, patchEdgeFaceRegion&);
+};
+
+
+//- Data associated with patchEdgeFaceRegion type are contiguous
+template<>
+inline bool contiguous<patchEdgeFaceRegion>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "patchEdgeFaceRegionI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceSources/regionToFace/patchEdgeFaceRegionI.H b/src/meshTools/sets/faceSources/regionToFace/patchEdgeFaceRegionI.H
new file mode 100644
index 0000000000..60a6d13eb3
--- /dev/null
+++ b/src/meshTools/sets/faceSources/regionToFace/patchEdgeFaceRegionI.H
@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update this with w2 if w2 nearer to pt.
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceRegion::update
+(
+    const patchEdgeFaceRegion& w2,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (!w2.valid(td))
+    {
+        FatalErrorIn("patchEdgeFaceRegion::update(..)")
+            << "problem." << abort(FatalError);
+    }
+
+    if (w2.region_ == -2 || region_ == -2)
+    {
+
+Pout<< "update : " << *this << "   w2:" << w2 << " return FALSE" << endl;
+
+        // Blocked edge/face
+        return false;
+    }
+
+    if (!valid(td))
+    {
+        // current not yet set so use any value
+        label oldRegion = region_;
+        operator=(w2);
+Pout<< "update : " << *this << " was:" << oldRegion
+    << "   w2:" << w2 << " return TRUE" << endl;
+        return true;
+    }
+    else
+    {
+        if (w2.region_ < region_)
+        {
+            label oldRegion = region_;
+            operator=(w2);
+Pout<< "update : " << *this << " was:" << oldRegion
+    << "   w2:" << w2 << " return TRUE" << endl;
+        return true;
+            return true;
+        }
+        else
+        {
+Pout<< "update : " << *this
+    << "   w2:" << w2 << " return FALSE" << endl;
+            return false;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::patchEdgeFaceRegion::patchEdgeFaceRegion()
+:
+    region_(-1)
+{}
+
+
+// Construct from origin, distance
+inline Foam::patchEdgeFaceRegion::patchEdgeFaceRegion
+(
+    const label region
+)
+:
+    region_(region)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::patchEdgeFaceRegion::region() const
+{
+    return region_;
+}
+
+
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceRegion::valid(TrackingData& td) const
+{
+    return region_ != -1;
+}
+
+
+template<class TrackingData>
+inline void Foam::patchEdgeFaceRegion::transform
+(
+    const polyMesh& mesh,
+    const indirectPrimitivePatch& patch,
+    const tensor& rotTensor,
+    const scalar tol,
+    TrackingData& td
+)
+{}
+
+
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceRegion::updateEdge
+(
+    const polyMesh& mesh,
+    const indirectPrimitivePatch& patch,
+    const label edgeI,
+    const label faceI,
+    const patchEdgeFaceRegion& faceInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    return update(faceInfo, tol, td);
+}
+
+
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceRegion::updateEdge
+(
+    const polyMesh& mesh,
+    const indirectPrimitivePatch& patch,
+    const patchEdgeFaceRegion& edgeInfo,
+    const bool sameOrientation,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    return update(edgeInfo, tol, td);
+}
+
+
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceRegion::updateFace
+(
+    const polyMesh& mesh,
+    const indirectPrimitivePatch& patch,
+    const label faceI,
+    const label edgeI,
+    const patchEdgeFaceRegion& edgeInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    return update(edgeInfo, tol, td);
+}
+
+
+template <class TrackingData>
+inline bool Foam::patchEdgeFaceRegion::equal
+(
+    const patchEdgeFaceRegion& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::patchEdgeFaceRegion::operator==
+(
+    const Foam::patchEdgeFaceRegion& rhs
+) const
+{
+    return region() == rhs.region();
+}
+
+
+inline bool Foam::patchEdgeFaceRegion::operator!=
+(
+    const Foam::patchEdgeFaceRegion& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceSources/regionToFace/regionToFace.C b/src/meshTools/sets/faceSources/regionToFace/regionToFace.C
new file mode 100644
index 0000000000..5c5f538876
--- /dev/null
+++ b/src/meshTools/sets/faceSources/regionToFace/regionToFace.C
@@ -0,0 +1,254 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionToFace.H"
+#include "polyMesh.H"
+#include "faceSet.H"
+#include "mappedPatchBase.H"
+#include "indirectPrimitivePatch.H"
+#include "PatchTools.H"
+#include "addToRunTimeSelectionTable.H"
+#include "PatchEdgeFaceWave.H"
+#include "patchEdgeFaceRegion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(regionToFace, 0);
+
+addToRunTimeSelectionTable(topoSetSource, regionToFace, word);
+
+addToRunTimeSelectionTable(topoSetSource, regionToFace, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::regionToFace::usage_
+(
+    regionToFace::typeName,
+    "\n    Usage: regionToFace <faceSet> (x y z)\n\n"
+    "    Select all faces in the connected region of the faceSet"
+    " starting from the point.\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::regionToFace::markZone
+(
+    const indirectPrimitivePatch& patch,
+    const label procI,
+    const label faceI,
+    const label zoneI,
+    labelList& faceZone
+) const
+{
+    // Data on all edges and faces
+    List<patchEdgeFaceRegion> allEdgeInfo(patch.nEdges());
+    List<patchEdgeFaceRegion> allFaceInfo(patch.size());
+
+    DynamicList<label> changedEdges;
+    DynamicList<patchEdgeFaceRegion> changedInfo;
+
+    if (Pstream::myProcNo() == procI)
+    {
+        const labelList& fEdges = patch.faceEdges()[faceI];
+        forAll(fEdges, i)
+        {
+            changedEdges.append(fEdges[i]);
+            changedInfo.append(zoneI);
+        }
+    }
+
+    // Walk
+    PatchEdgeFaceWave
+    <
+        indirectPrimitivePatch,
+        patchEdgeFaceRegion
+    > calc
+    (
+        mesh_,
+        patch,
+        changedEdges,
+        changedInfo,
+        allEdgeInfo,
+        allFaceInfo,
+        returnReduce(patch.nEdges(), sumOp<label>())
+    );
+
+    forAll(allFaceInfo, faceI)
+    {
+        if (allFaceInfo[faceI].region() == zoneI)
+        {
+            faceZone[faceI] = zoneI;
+        }
+    }
+}
+
+
+void Foam::regionToFace::combine(topoSet& set, const bool add) const
+{
+    Info<< "    Loading subset " << setName_ << " to delimit search region."
+        << endl;
+    faceSet subSet(mesh_, setName_);
+
+    indirectPrimitivePatch patch
+    (
+        IndirectList<face>(mesh_.faces(), subSet.toc()),
+        mesh_.points()
+    );
+
+    mappedPatchBase::nearInfo ni
+    (
+        pointIndexHit(false, vector::zero, -1),
+        Tuple2<scalar, label>
+        (
+            sqr(GREAT),
+            Pstream::myProcNo()
+        )
+    );
+
+    forAll(patch, i)
+    {
+        const point& fc = patch.faceCentres()[i];
+        scalar d2 = magSqr(fc-nearPoint_);
+
+        if (!ni.first().hit() || d2 < ni.second().first())
+        {
+            ni.second().first() = d2;
+            ni.first().setHit();
+            ni.first().setPoint(fc);
+            ni.first().setIndex(i);
+        }
+    }
+
+    // Globally reduce
+    combineReduce(ni, mappedPatchBase::nearestEqOp());
+
+    Info<< "    Found nearest face at " << ni.first().rawPoint()
+        << " on processor " << ni.second().second()
+        << " face " << ni.first().index()
+        << " distance " << Foam::sqrt(ni.second().first()) << endl;
+
+    labelList faceRegion(patch.size(), -1);
+    markZone
+    (
+        patch,
+        ni.second().second(),   // procI
+        ni.first().index(),     // start face
+        0,                      // currentZone
+        faceRegion
+    );
+
+    forAll(faceRegion, faceI)
+    {
+        if (faceRegion[faceI] == 0)
+        {
+            addOrDelete(set, patch.addressing()[faceI], add);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::regionToFace::regionToFace
+(
+    const polyMesh& mesh,
+    const word& setName,
+    const point& nearPoint
+)
+:
+    topoSetSource(mesh),
+    setName_(setName),
+    nearPoint_(nearPoint)
+{}
+
+
+// Construct from dictionary
+Foam::regionToFace::regionToFace
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    setName_(dict.lookup("set")),
+    nearPoint_(dict.lookup("nearPoint"))
+{}
+
+
+// Construct from Istream
+Foam::regionToFace::regionToFace
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    setName_(checkIs(is)),
+    nearPoint_(checkIs(is))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::regionToFace::~regionToFace()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::regionToFace::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+    {
+        Info<< "    Adding all faces of connected region of set "
+            << setName_
+            << " starting from point "
+            << nearPoint_ << " ..." << endl;
+
+        combine(set, true);
+    }
+    else if (action == topoSetSource::DELETE)
+    {
+        Info<< "    Removing all cells of connected region of set "
+            << setName_
+            << " starting from point "
+            << nearPoint_ << " ..." << endl;
+
+        combine(set, false);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceSources/regionToFace/regionToFace.H b/src/meshTools/sets/faceSources/regionToFace/regionToFace.H
new file mode 100644
index 0000000000..1dc1d48069
--- /dev/null
+++ b/src/meshTools/sets/faceSources/regionToFace/regionToFace.H
@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionToFace
+
+Description
+    A topoSetSource to select faces belonging to topological connected region
+    (that contains given point)
+
+SourceFiles
+    regionToFace.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionToFace_H
+#define regionToFace_H
+
+#include "topoSetSource.H"
+#include "PackedBoolList.H"
+#include "indirectPrimitivePatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class regionToFace Declaration
+\*---------------------------------------------------------------------------*/
+
+class regionToFace
+:
+    public topoSetSource
+{
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name of set to use
+        word setName_;
+
+        //- Coordinate that is nearest/on connected region
+        point nearPoint_;
+
+    // Private Member Functions
+
+        //- Walk edge-face-edge
+        void markZone
+        (
+            const indirectPrimitivePatch& patch,
+            const label procI,
+            const label faceI,
+            const label zoneI,
+            labelList& faceZone
+        ) const;
+
+        void combine(topoSet& set, const bool add) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("regionToFace");
+
+    // Constructors
+
+        //- Construct from components
+        regionToFace
+        (
+            const polyMesh& mesh,
+            const word& setName,
+            const point& nearPoint
+        );
+
+        //- Construct from dictionary
+        regionToFace
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        regionToFace
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    //- Destructor
+    virtual ~regionToFace();
+
+
+    // Member Functions
+
+        virtual sourceType setType() const
+        {
+            return FACESETSOURCE;
+        }
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
index ecedee0485..1cb392eb6a 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
@@ -255,6 +255,9 @@ void thermoSingleLayer::updateSubmodels()
     // Update source fields
     hsSp_ += primaryEnergyPCTrans_/magSf()/time().deltaT();
     rhoSp_ += primaryMassPCTrans_/magSf()/time().deltaT();
+
+    // Vapour recoil pressure
+    pSp_ -= sqr(primaryMassPCTrans_/magSf()/time_.deltaT())/2.0/rhoPrimary_;
 }
 
 
diff --git a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/combustionProperties b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/combustionProperties
index 41c94a4fc1..6c1f0e9cf3 100644
--- a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/combustionProperties
+++ b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/combustionProperties
@@ -15,18 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<psiChemistryCombustionModel>;
+combustionModel  PaSR<psiChemistryCombustion>;
 
 active  true;
 
-infinitelyFastChemistryCoeffs
-{
-    C       10.0;
-}
-
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 0.1;
+    Cmix                0.1;
     turbulentReaction   on;
 }
 
diff --git a/tutorials/incompressible/icoFoam/elbow/system/foamDataToFluentDict b/tutorials/incompressible/icoFoam/elbow/system/foamDataToFluentDict
index 2e4c170627..56e694db05 100644
--- a/tutorials/incompressible/icoFoam/elbow/system/foamDataToFluentDict
+++ b/tutorials/incompressible/icoFoam/elbow/system/foamDataToFluentDict
@@ -27,7 +27,7 @@ k               5;
 
 epsilon         6;
 
-gamma           150;
+alpha1          150;
 
 
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/constant/combustionProperties b/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/constant/combustionProperties
index e531cc4a16..754f22dad2 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/constant/combustionProperties
+++ b/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<rhoChemistryCombustionModel>;
+combustionModel  PaSR<rhoChemistryCombustion>;
 
 active  true;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 0.1;
+    Cmix                0.1;
     turbulentReaction   on;
     useReactionRate     true;
 }
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/combustionProperties b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/combustionProperties
index 6ab266d91f..5ede7572c3 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/combustionProperties
+++ b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<rhoChemistryCombustionModel>;
+combustionModel  PaSR<rhoChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 0.1;
+    Cmix                0.1;
     turbulentReaction   off;
     useReactionRate     true;
 }
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/combustionProperties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/combustionProperties
index d2fe2410f0..1832dddf0a 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/combustionProperties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/combustionProperties
@@ -15,14 +15,14 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<psiChemistryCombustionModel>;
+combustionModel  PaSR<psiChemistryCombustion>;
 
 active  true;
 
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   on;
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/combustionProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/combustionProperties
index cec49a34e9..334391590c 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/combustionProperties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<rhoChemistryCombustionModel>;
+combustionModel  PaSR<rhoChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   off;
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/combustionProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/combustionProperties
index cec49a34e9..334391590c 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/combustionProperties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<rhoChemistryCombustionModel>;
+combustionModel  PaSR<rhoChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   off;
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/combustionProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/combustionProperties
index cec49a34e9..334391590c 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/combustionProperties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<rhoChemistryCombustionModel>;
+combustionModel  PaSR<rhoChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   off;
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/combustionProperties b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/combustionProperties
index acda0a53d4..4fe87a3132 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/combustionProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<psiChemistryCombustionModel>;
+combustionModel  PaSR<psiChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   on;
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/combustionProperties b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/combustionProperties
index 3a7e78cb55..87d04695b8 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/combustionProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<psiChemistryCombustionModel>;
+combustionModel  PaSR<psiChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   off;
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/combustionProperties b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/combustionProperties
index 3a7e78cb55..87d04695b8 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/combustionProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<psiChemistryCombustionModel>;
+combustionModel  PaSR<psiChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   off;
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/combustionProperties b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/combustionProperties
index 3a7e78cb55..87d04695b8 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/combustionProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<psiChemistryCombustionModel>;
+combustionModel  PaSR<psiChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   off;
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/combustionProperties b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/combustionProperties
index 3a7e78cb55..87d04695b8 100644
--- a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/combustionProperties
+++ b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<psiChemistryCombustionModel>;
+combustionModel  PaSR<psiChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   off;
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes
index dc13503bc6..a59829fcd8 100644
--- a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes
+++ b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes
@@ -32,12 +32,11 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
     div(phi,K)      Gauss linear;
-    div(phi,h)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(U)          Gauss linear;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
-    div(phi,Yi_h)   Gauss upwind;
+    div(phi,Yi_hs)  Gauss upwind;
 }
 
 laplacianSchemes
@@ -49,7 +48,6 @@ laplacianSchemes
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(DREff,R) Gauss linear corrected;
     laplacian((rho*(1|A(U))),p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
 }
 
 interpolationSchemes
diff --git a/tutorials/lagrangian/sprayFoam/aachenBomb/constant/combustionProperties b/tutorials/lagrangian/sprayFoam/aachenBomb/constant/combustionProperties
index 906a38d81b..8f326f3667 100644
--- a/tutorials/lagrangian/sprayFoam/aachenBomb/constant/combustionProperties
+++ b/tutorials/lagrangian/sprayFoam/aachenBomb/constant/combustionProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  PaSR<psiChemistryCombustionModel>;
+combustionModel  PaSR<psiChemistryCombustion>;
 
 active  false;
 
 PaSRCoeffs
 {
-    Cmix                Cmix [ 0 0 0 0 0 0 0 ] 1.0;
+    Cmix                1.0;
     turbulentReaction   yes;
 }
 
diff --git a/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes b/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes
index 74c85361d1..c9c18972f1 100644
--- a/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes
+++ b/tutorials/lagrangian/sprayFoam/aachenBomb/system/fvSchemes
@@ -32,12 +32,11 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
     div(phi,K)      Gauss linear;
-    div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(U)          Gauss linear;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
-    div(phi,Yi_h)   Gauss upwind;
+    div(phi,Yi_hs)  Gauss upwind;
 }
 
 laplacianSchemes