diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
index 365e64d1b1..ca2936633a 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
@@ -1,5 +1,8 @@
 EXE_INC = \
     -I../buoyantBoussinesqSimpleFoam \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
@@ -9,6 +12,8 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfvOptions \
+    -lsampling \
     -lmeshTools \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H
index 65c92dab51..76b0ef0168 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H
@@ -11,12 +11,18 @@
       - fvm::laplacian(alphaEff, T)
      ==
         radiation->ST(rhoCpRef, T)
+      + fvOptions(T)
     );
 
     TEqn.relax();
+
+    fvOptions.constrain(TEqn);
+
     TEqn.solve();
 
     radiation->correct();
 
+    fvOptions.correct(T);
+
     rhok = 1.0 - beta*(T - TRef);
 }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
index 07f46ec998..92768088b9 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
@@ -5,10 +5,14 @@
         fvm::ddt(U)
       + fvm::div(phi, U)
       + turbulence->divDevReff(U)
+     ==
+        fvOptions(U)
     );
 
     UEqn.relax();
 
+    fvOptions.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve
@@ -23,4 +27,6 @@
                 )*mesh.magSf()
             )
         );
+
+        fvOptions.correct(U);
     }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index 65dda0653e..6168082e1d 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -48,8 +48,9 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
-#include "pimpleControl.H"
 #include "radiationModel.H"
+#include "fvIOoptionList.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,6 +62,7 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createIncompressibleRadiationModel.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "CourantNo.H"
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options
index 6905967f04..6a6311e47c 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options
@@ -1,5 +1,8 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
     -I$(LIB_SRC)/transportModels \
@@ -7,7 +10,9 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lsampling \
     -lmeshTools \
+    -lfvOptions \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleTransportModels
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
index a0c5d24b2b..c495e58285 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
@@ -8,10 +8,17 @@
     (
         fvm::div(phi, T)
       - fvm::laplacian(alphaEff, T)
+     ==
+        fvOptions(T)
     );
 
     TEqn.relax();
+
+    fvOptions.constrain(TEqn);
+
     TEqn.solve();
 
+    fvOptions.correct(T);
+
     rhok = 1.0 - beta*(T - TRef);
 }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
index cbe464fc02..dd516edf2f 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
@@ -4,16 +4,20 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
+     ==
+        fvOptions(U)
     );
 
     UEqn().relax();
 
+    fvOptions.constrain(UEqn());
+
     if (simple.momentumPredictor())
     {
         solve
         (
             UEqn()
-            ==
+          ==
             fvc::reconstruct
             (
                 (
@@ -22,4 +26,6 @@
                 )*mesh.magSf()
             )
         );
+
+        fvOptions.correct(U);
     }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
index 4fc37efd28..b0d2f98aeb 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
@@ -48,6 +48,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "fvIOoptionList.H"
 #include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -59,6 +60,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     simpleControl simple(mesh);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index 5297be327f..b2cd9ac54e 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -26,21 +26,18 @@ Application
 
 Description
     Combination of heatConductionFoam and buoyantFoam for conjugate heat
-    transfer between a solid region and fluid region. It includes
-    porous media in the primary fluid region treated explicitly.
+    transfer between solid regions and fluid regions. Both regions include
+    the fvOptions framework.
 
     It handles secondary fluid or solid circuits which can be coupled
     thermally with the main fluid region. i.e radiators, etc.
 
-    The secondary fluid region is
-
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "rhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
-#include "zeroGradientFvPatchFields.H"
 #include "regionProperties.H"
 #include "compressibleCourantNo.H"
 #include "solidRegionDiffNo.H"
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
index 000948514e..de963fe3b2 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
@@ -30,7 +30,7 @@ Description
 
     Sub-models include:
     - turbulence modelling, i.e. laminar, RAS or LES
-    - run-time selectable finitie volume options, e.g. MRF, explicit porosity
+    - run-time selectable finite volume options, e.g. MRF, explicit porosity
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/applications/utilities/postProcessing/velocityField/Pe/Pe.C b/applications/utilities/postProcessing/velocityField/Pe/Pe.C
index 0c4f61729c..544a116921 100644
--- a/applications/utilities/postProcessing/velocityField/Pe/Pe.C
+++ b/applications/utilities/postProcessing/velocityField/Pe/Pe.C
@@ -130,7 +130,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                             mesh.magSf()
                           * mesh.surfaceInterpolation::deltaCoeffs()
                           * fvc::interpolate(RASModel->nuEff())
-                         )
+                        )
                     )
                 );
             }
@@ -192,8 +192,12 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                             mesh,
                             IOobject::NO_READ
                         ),
-                        mesh.surfaceInterpolation::deltaCoeffs()
-                      * (mag(phi)/mesh.magSf())*(runTime.deltaT()/nu)
+                        mag(phi)
+                       /(
+                            mesh.magSf()
+                          * mesh.surfaceInterpolation::deltaCoeffs()
+                          * nu
+                        )
                     )
                 );
             }
@@ -317,8 +321,12 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                             mesh,
                             IOobject::NO_READ
                         ),
-                        mesh.surfaceInterpolation::deltaCoeffs()
-                      * (mag(phi)/(mesh.magSf()))*(runTime.deltaT()/mu)
+                        mag(phi)
+                       /(
+                            mesh.magSf()
+                          * mesh.surfaceInterpolation::deltaCoeffs()
+                          * mu
+                        )
                     )
                 );
             }
@@ -330,26 +338,6 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                     << abort(FatalError);
         }
 
-
-        // can also check how many cells exceed a particular Pe limit
-        /*
-        {
-            label count = 0;
-            label PeLimit = 200;
-            forAll(PePtr(), i)
-            {
-                if (PePtr()[i] > PeLimit)
-                {
-                    count++;
-                }
-
-            }
-
-            Info<< "Fraction > " << PeLimit << " = "
-                << scalar(count)/Pe.size() << endl;
-        }
-        */
-
         Info<< "Pe max : " << max(PePtr()).value() << endl;
 
         if (writeResults)
diff --git a/applications/utilities/surface/surfaceCheck/surfaceCheck.C b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
index a017055958..29707e303c 100644
--- a/applications/utilities/surface/surfaceCheck/surfaceCheck.C
+++ b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
@@ -653,6 +653,36 @@ int main(int argc, char *argv[])
         Info<< "Checking self-intersection." << endl;
 
         triSurfaceSearch querySurf(surf);
+
+        //{
+        //    OBJstream intStream("selfInter2.obj");
+        //    const indexedOctree<treeDataTriSurface>& tree = querySurf.tree();
+        //    forAll(surf.edges(), edgeI)
+        //    {
+        //        const edge& e = surf.edges()[edgeI];
+        //
+        //        pointIndexHit hitInfo
+        //        (
+        //            tree.findLine
+        //            (
+        //                surf.points()[surf.meshPoints()[e[0]]],
+        //                surf.points()[surf.meshPoints()[e[1]]],
+        //                treeDataTriSurface::findSelfIntersectOp
+        //                (
+        //                    tree,
+        //                    edgeI
+        //                )
+        //            )
+        //        );
+        //
+        //        if (hitInfo.hit())
+        //        {
+        //            Pout<< "Found hit:" << hitInfo.hitPoint() << endl;
+        //            intStream.write(hitInfo.hitPoint());
+        //        }
+        //    }
+        //}
+
         surfaceIntersection inter(querySurf);
 
         if (inter.cutEdges().empty() && inter.cutPoints().empty())
diff --git a/etc/codeTemplates/dynamicCode/functionObjectTemplate.C b/etc/codeTemplates/dynamicCode/functionObjectTemplate.C
index b44bdfb8dc..91b31dda2b 100644
--- a/etc/codeTemplates/dynamicCode/functionObjectTemplate.C
+++ b/etc/codeTemplates/dynamicCode/functionObjectTemplate.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -123,6 +123,19 @@ void ${typeName}FunctionObject::end()
 }
 
 
+void ${typeName}FunctionObject::timeSet()
+{
+    if (${verbose:-false})
+    {
+        Info<<"timeSet ${typeName} sha1: ${SHA1sum}\n";
+    }
+
+//{{{ begin codeTime
+    ${codeTimeSet}
+//}}} end code
+}
+
+
 void ${typeName}FunctionObject::write()
 {
     if (${verbose:-false})
diff --git a/etc/codeTemplates/dynamicCode/functionObjectTemplate.H b/etc/codeTemplates/dynamicCode/functionObjectTemplate.H
index 211c0868ec..a139f25a8b 100644
--- a/etc/codeTemplates/dynamicCode/functionObjectTemplate.H
+++ b/etc/codeTemplates/dynamicCode/functionObjectTemplate.H
@@ -127,6 +127,9 @@ public:
         //- Write, execute the "writeCalls"
         virtual void write();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Update for changes of mesh
         virtual void updateMesh(const mapPolyMesh&)
         {}
diff --git a/src/OpenFOAM/fields/Fields/Field/Field.C b/src/OpenFOAM/fields/Fields/Field/Field.C
index 334eba8e41..7460dba215 100644
--- a/src/OpenFOAM/fields/Fields/Field/Field.C
+++ b/src/OpenFOAM/fields/Fields/Field/Field.C
@@ -677,6 +677,7 @@ template<class Type>
 template<class Form, class Cmpt, int nCmpt>
 void Foam::Field<Type>::operator=(const VectorSpace<Form,Cmpt,nCmpt>& vs)
 {
+    typedef VectorSpace<Form,Cmpt,nCmpt> VSType;
     TFOR_ALL_F_OP_S(Type, *this, =, VSType, vs)
 }
 
diff --git a/src/OpenFOAM/fields/Fields/scalarField/scalarField.C b/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
index 96e41d9571..7ae37531de 100644
--- a/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
+++ b/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "scalarField.H"
+#include "unitConversion.H"
 
 #define TEMPLATE
 #include "FieldFunctionsM.C"
@@ -155,6 +156,11 @@ UNARY_FUNCTION(scalar, scalar, j1)
 UNARY_FUNCTION(scalar, scalar, y0)
 UNARY_FUNCTION(scalar, scalar, y1)
 
+UNARY_FUNCTION(scalar, scalar, degToRad)
+UNARY_FUNCTION(scalar, scalar, radToDeg)
+UNARY_FUNCTION(scalar, scalar, atmToPa)
+UNARY_FUNCTION(scalar, scalar, paToAtm)
+
 
 #define BesselFunc(func)                                                      \
 void func(scalarField& res, const int n, const UList<scalar>& sf)             \
diff --git a/src/OpenFOAM/fields/Fields/scalarField/scalarField.H b/src/OpenFOAM/fields/Fields/scalarField/scalarField.H
index ff322d3fe4..f30938e9b6 100644
--- a/src/OpenFOAM/fields/Fields/scalarField/scalarField.H
+++ b/src/OpenFOAM/fields/Fields/scalarField/scalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -130,6 +130,10 @@ UNARY_FUNCTION(scalar, scalar, j1)
 UNARY_FUNCTION(scalar, scalar, y0)
 UNARY_FUNCTION(scalar, scalar, y1)
 
+UNARY_FUNCTION(scalar, scalar, degToRad)
+UNARY_FUNCTION(scalar, scalar, radToDeg)
+UNARY_FUNCTION(scalar, scalar, atmToPa)
+UNARY_FUNCTION(scalar, scalar, paToAtm)
 
 #define BesselFunc(func)                                            \
 void func(scalarField& Res, const int n, const UList<scalar>& sf);  \
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
index 7f23facd04..ef500b4b8b 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
@@ -63,16 +63,6 @@ public:
         const vectorField& cellCtrs
     );
 
-//    static tmp<scalarField> faceOrthogonality(const polyMesh& mesh)
-//    {
-//        return faceOrthogonality
-//        (
-//            mesh,
-//            mesh.faceAreas(),
-//            mesh.cellCentres()
-//        );
-//    }
-
     //- Generate skewness field
     static tmp<scalarField> faceSkewness
     (
@@ -83,18 +73,6 @@ public:
         const vectorField& cellCtrs
     );
 
-//    static tmp<scalarField> faceSkewness(const polyMesh& mesh)
-//    {
-//        return faceSkewness
-//        (
-//            mesh,
-//            mesh.points(),
-//            mesh.faceCentres(),
-//            mesh.faceAreas(),
-//            mesh.cellCentres()
-//        );
-//    }
-
     //- Generate interpolation factors field
     static tmp<scalarField> faceWeights
     (
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
index 738f87d0a8..81d722a57b 100644
--- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
@@ -88,7 +88,17 @@ Foam::porosityModel::porosityModel
     const word& cellZoneName
 )
 :
-    MeshObject<fvMesh, Foam::UpdateableMeshObject, porosityModel>(mesh),
+    regIOobject
+    (
+        IOobject
+        (
+            name,
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    ),
     name_(name),
     mesh_(mesh),
     dict_(dict),
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H
index 36a5a9eafb..06df0bfdae 100644
--- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H
@@ -36,7 +36,6 @@ SourceFiles
 #ifndef porosityModel_H
 #define porosityModel_H
 
-#include "MeshObject.H"
 #include "fvMesh.H"
 #include "dictionary.H"
 #include "fvMatricesFwd.H"
@@ -56,7 +55,7 @@ namespace Foam
 
 class porosityModel
 :
-    public MeshObject<fvMesh, UpdateableMeshObject, porosityModel>
+    public regIOobject
 {
 private:
 
diff --git a/src/meshTools/indexedOctree/treeDataEdge.C b/src/meshTools/indexedOctree/treeDataEdge.C
index 4d0b9b4d81..2499380b32 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.C
+++ b/src/meshTools/indexedOctree/treeDataEdge.C
@@ -109,8 +109,6 @@ Foam::treeDataEdge::findIntersectOp::findIntersectOp
 (
     const indexedOctree<treeDataEdge>& tree
 )
-:
-    tree_(tree)
 {}
 
 
diff --git a/src/meshTools/indexedOctree/treeDataEdge.H b/src/meshTools/indexedOctree/treeDataEdge.H
index 12fd0d6aa1..be82a2a09c 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.H
+++ b/src/meshTools/indexedOctree/treeDataEdge.H
@@ -121,8 +121,6 @@ public:
 
     class findIntersectOp
     {
-        const indexedOctree<treeDataEdge>& tree_;
-
     public:
 
         findIntersectOp(const indexedOctree<treeDataEdge>& tree);
diff --git a/src/meshTools/indexedOctree/treeDataPoint.C b/src/meshTools/indexedOctree/treeDataPoint.C
index fd7f4404b4..82a9bd257b 100644
--- a/src/meshTools/indexedOctree/treeDataPoint.C
+++ b/src/meshTools/indexedOctree/treeDataPoint.C
@@ -70,8 +70,6 @@ Foam::treeDataPoint::findIntersectOp::findIntersectOp
 (
     const indexedOctree<treeDataPoint>& tree
 )
-:
-    tree_(tree)
 {}
 
 
diff --git a/src/meshTools/indexedOctree/treeDataPoint.H b/src/meshTools/indexedOctree/treeDataPoint.H
index 020f60a11d..73cec8a212 100644
--- a/src/meshTools/indexedOctree/treeDataPoint.H
+++ b/src/meshTools/indexedOctree/treeDataPoint.H
@@ -104,8 +104,6 @@ public:
 
     class findIntersectOp
     {
-        const indexedOctree<treeDataPoint>& tree_;
-
     public:
 
         findIntersectOp(const indexedOctree<treeDataPoint>& tree);
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
index 488104a345..51e5252121 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
@@ -66,76 +66,6 @@ void Foam::treeDataPrimitivePatch<PatchType>::update()
 }
 
 
-template<class PatchType>
-bool Foam::treeDataPrimitivePatch<PatchType>::findIntersection
-(
-    const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
-    const label index,
-    const point& start,
-    const point& end,
-    point& intersectionPoint
-)
-{
-    const treeDataPrimitivePatch<PatchType>& shape = tree.shapes();
-    const PatchType& patch = shape.patch();
-
-    const pointField& points = patch.points();
-    const typename PatchType::FaceType& f = patch[index];
-
-    // Do quick rejection test
-    if (shape.cacheBb_)
-    {
-        const treeBoundBox& faceBb = shape.bbs_[index];
-
-        if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
-        {
-            // start and end in same block outside of faceBb.
-            return false;
-        }
-    }
-
-    const vector dir(end - start);
-    pointHit inter;
-
-    if (f.size() == 3)
-    {
-        inter = triPointRef
-        (
-            points[f[0]],
-            points[f[1]],
-            points[f[2]]
-        ).intersection(start, dir, intersection::HALF_RAY, shape.planarTol_);
-    }
-    else
-    {
-        const pointField& faceCentres = patch.faceCentres();
-
-        inter = f.intersection
-        (
-            start,
-            dir,
-            faceCentres[index],
-            points,
-            intersection::HALF_RAY,
-            shape.planarTol_
-        );
-    }
-
-    if (inter.hit() && inter.distance() <= 1)
-    {
-        // Note: no extra test on whether intersection is in front of us
-        // since using half_ray
-        intersectionPoint = inter.hitPoint();
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
@@ -187,6 +117,19 @@ Foam::treeDataPrimitivePatch<PatchType>::findAllIntersectOp::findAllIntersectOp
 {}
 
 
+template<class PatchType>
+Foam::treeDataPrimitivePatch<PatchType>::
+findSelfIntersectOp::findSelfIntersectOp
+(
+    const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
+    const label edgeID
+)
+:
+    tree_(tree),
+    edgeID_(edgeID)
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class PatchType>
@@ -645,4 +588,116 @@ bool Foam::treeDataPrimitivePatch<PatchType>::findAllIntersectOp::operator()
 }
 
 
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::findSelfIntersectOp::operator()
+(
+    const label index,
+    const point& start,
+    const point& end,
+    point& intersectionPoint
+) const
+{
+    if (edgeID_ == -1)
+    {
+        FatalErrorIn
+        (
+            "findSelfIntersectOp::operator()\n"
+            "(\n"
+            "    const label index,\n"
+            "    const point& start,\n"
+            "    const point& end,\n"
+            "    point& intersectionPoint\n"
+            ") const"
+        )   << "EdgeID not set. Please set edgeID to the index of"
+            << " the edge you are testing"
+            << exit(FatalError);
+    }
+
+    const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
+    const PatchType& patch = shape.patch();
+
+    const typename PatchType::FaceType& f = patch.localFaces()[index];
+    const edge& e = patch.edges()[edgeID_];
+
+    if (findIndex(f, e[0]) == -1 && findIndex(f, e[1]) == -1)
+    {
+        return findIntersection(tree_, index, start, end, intersectionPoint);
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::findIntersection
+(
+    const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
+    const label index,
+    const point& start,
+    const point& end,
+    point& intersectionPoint
+)
+{
+    const treeDataPrimitivePatch<PatchType>& shape = tree.shapes();
+    const PatchType& patch = shape.patch();
+
+    const pointField& points = patch.points();
+    const typename PatchType::FaceType& f = patch[index];
+
+    // Do quick rejection test
+    if (shape.cacheBb_)
+    {
+        const treeBoundBox& faceBb = shape.bbs_[index];
+
+        if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
+        {
+            // start and end in same block outside of faceBb.
+            return false;
+        }
+    }
+
+    const vector dir(end - start);
+    pointHit inter;
+
+    if (f.size() == 3)
+    {
+        inter = triPointRef
+        (
+            points[f[0]],
+            points[f[1]],
+            points[f[2]]
+        ).intersection(start, dir, intersection::HALF_RAY, shape.planarTol_);
+    }
+    else
+    {
+        const pointField& faceCentres = patch.faceCentres();
+
+        inter = f.intersection
+        (
+            start,
+            dir,
+            faceCentres[index],
+            points,
+            intersection::HALF_RAY,
+            shape.planarTol_
+        );
+    }
+
+    if (inter.hit() && inter.distance() <= 1)
+    {
+        // Note: no extra test on whether intersection is in front of us
+        // since using half_ray
+        intersectionPoint = inter.hitPoint();
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
index 9d43d2ba85..18c901244e 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
@@ -86,16 +86,6 @@ class treeDataPrimitivePatch
         //- Initialise all member data
         void update();
 
-        //- Find intersection of line with shapes
-        static bool findIntersection
-        (
-            const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
-            const label index,
-            const point& start,
-            const point& end,
-            point& intersectionPoint
-        );
-
 
 public:
 
@@ -141,8 +131,8 @@ public:
 
         findIntersectOp(const indexedOctree<treeDataPrimitivePatch>& tree);
 
-        //- Calculate intersection of triangle with ray. Sets result
-        //  accordingly
+        //- Calculate intersection of any face with ray. Sets result
+        //  accordingly. Used to find first intersection.
         bool operator()
         (
             const label index,
@@ -167,8 +157,34 @@ public:
             DynamicList<label>& shapeMask
         );
 
-        //- Calculate intersection of triangle with ray. Sets result
-        //  accordingly
+        //- Calculate intersection of unique face with ray. Sets result
+        //  accordingly. Used to find all faces.
+        bool operator()
+        (
+            const label index,
+            const point& start,
+            const point& end,
+            point& intersectionPoint
+        ) const;
+    };
+
+
+    class findSelfIntersectOp
+    {
+        const indexedOctree<treeDataPrimitivePatch>& tree_;
+
+        const label edgeID_;
+
+    public:
+
+        findSelfIntersectOp
+        (
+            const indexedOctree<treeDataPrimitivePatch>& tree,
+            const label edgeID
+        );
+
+        //- Calculate intersection of face with edge of patch. Excludes
+        //  faces that use edgeID. Used to find self intersection.
         bool operator()
         (
             const label index,
@@ -234,6 +250,16 @@ public:
                 const point& centre,
                 const scalar radiusSqr
             ) const;
+
+            //- Helper: find intersection of line with shapes
+            static bool findIntersection
+            (
+                const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
+                const label index,
+                const point& start,
+                const point& end,
+                point& intersectionPoint
+            );
 };
 
 
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 7c80bd1a56..95d0e887b6 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -74,6 +74,26 @@ const Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>
 // Read my additional data from the dictionary
 bool Foam::distributedTriSurfaceMesh::read()
 {
+
+    if
+    (
+        Pstream::parRun()
+     &&
+     (
+        regIOobject::fileModificationChecking == timeStampMaster
+     || regIOobject::fileModificationChecking == inotifyMaster
+     )
+    )
+    {
+        FatalErrorIn("Foam::distributedTriSurfaceMesh::read()")
+            << "    distributedTriSurfaceMesh is being constructed\n"
+            << "    using 'timeStampMaster' or 'inotifyMaster.'\n"
+            << "    Modify the entry fileModificationChecking\n"
+            << "    in the etc/controlDict.\n"
+            << "    Use 'timeStamp' instead."
+            << exit(FatalError);
+    }
+
     // Get bb of all domains.
     procBb_.setSize(Pstream::nProcs());
 
@@ -183,10 +203,6 @@ void Foam::distributedTriSurfaceMesh::distributeSegment
     List<DynamicList<label> >& sendMap
 ) const
 {
-    // Work points
-    point clipPt;
-
-
     // 1. Fully local already handled outside. Note: retest is cheap.
     if (isLocal(procBb_[Pstream::myProcNo()], start, end))
     {
@@ -227,6 +243,9 @@ void Foam::distributedTriSurfaceMesh::distributeSegment
             // Scheme a: any processor that intersects the segment gets
             // the segment.
 
+            // Intersection point
+            point clipPt;
+
             if (bb.intersects(start, end, clipPt))
             {
                 sendMap[procI].append(allSegments.size());
diff --git a/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C b/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C
index 5c1ed4a967..cd8c3e0169 100644
--- a/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C
+++ b/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C
@@ -61,6 +61,7 @@ void Foam::codedFunctionObject::prepare
     dynCode.setFilterVariable("codeExecute", codeExecute_);
     dynCode.setFilterVariable("codeEnd", codeEnd_);
     dynCode.setFilterVariable("codeData", codeData_);
+    dynCode.setFilterVariable("codeTimeSet", codeTimeSet_);
     //dynCode.setFilterVariable("codeWrite", codeWrite_);
 
     // compile filtered C template
@@ -187,7 +188,8 @@ bool Foam::codedFunctionObject::end()
 
 bool Foam::codedFunctionObject::timeSet()
 {
-    return false;
+    updateLibrary(redirectType_);
+    return redirectFunctionObject().timeSet();
 }
 
 
@@ -267,6 +269,24 @@ bool Foam::codedFunctionObject::read(const dictionary& dict)
         );
     }
 
+    const entry* timeSetPtr = dict.lookupEntryPtr
+    (
+        "codeTimeSet",
+        false,
+        false
+    );
+    if (timeSetPtr)
+    {
+        codeTimeSet_ = stringOps::trim(timeSetPtr->stream());
+        stringOps::inplaceExpand(codeTimeSet_, dict);
+        dynamicCodeContext::addLineDirective
+        (
+            codeTimeSet_,
+            timeSetPtr->startLineNumber(),
+            dict.name()
+        );
+    }
+
     updateLibrary(redirectType_);
     return redirectFunctionObject().read(dict);
 }
diff --git a/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.H b/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.H
index f597407e1c..7b3f8f57e3 100644
--- a/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.H
+++ b/src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.H
@@ -31,6 +31,39 @@ Description
     This function object provides a general interface to enable dynamic code
     compilation.
 
+    The entries are
+        code        : c++; upon functionObject::write()
+        codeInclude : include files
+        codeOptions : include paths; inserted into EXE_INC in Make/options
+        codeLibs    : link line; inserted into LIB_LIBS in Make/options
+
+        codeExecute : c++;upon functionObject::execute();
+        codeRead    : c++; upon functionObject::read();
+        codeEnd     : c++; upon functionObject::end();
+        codeData    : c++; local member data (null constructed);
+        codeTimeSet : c++; upon functionObject::timeSet();
+        localCode   : c++; local static functions
+
+    Example of function object specification:
+    \verbatim
+    difference
+    {
+        functionObjectLibs ("libutilityFunctionObjects.so");
+
+        type coded;
+        // Name of on-the-fly generated functionObject
+        redirectType writeMagU;
+        code
+        #{
+            // Lookup U
+            const volVectorField& U = mesh().lookupObject<volVectorField>("U");
+            // Write
+            mag(U).write();
+        }
+    }
+    \endverbatim
+
+
 SeeAlso
     Foam::functionObject
     Foam::OutputFilterFunctionObject
@@ -77,6 +110,7 @@ protected:
         string codeRead_;
         string codeExecute_;
         string codeEnd_;
+        string codeTimeSet_;
 
         //- Underlying functionObject
         mutable autoPtr<functionObject> redirectFunctionObjectPtr_;
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C
index 108109702b..bc4bf5ee88 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C
@@ -76,7 +76,7 @@ tmp<volVectorField> laminar::Us() const
         (
             IOobject
             (
-                typeName + ".Us",
+                typeName + ":Us",
                 owner_.regionMesh().time().timeName(),
                 owner_.regionMesh(),
                 IOobject::NO_READ,
@@ -104,7 +104,7 @@ tmp<volScalarField> laminar::mut() const
         (
             IOobject
             (
-                typeName + ".mut",
+                typeName + ":mut",
                 owner_.regionMesh().time().timeName(),
                 owner_.regionMesh(),
                 IOobject::NO_READ,
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
index 4a08755ee3..8b0a0c5f16 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,8 +26,6 @@ License
 #include "standardPhaseChange.H"
 #include "addToRunTimeSelectionTable.H"
 #include "thermoSingleLayer.H"
-#include "specie.H"
-#include "heatTransferModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
index 42d459ce01..03d1654fda 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
@@ -294,10 +294,20 @@ tmp<fvScalarMatrix> thermoSingleLayer::q(volScalarField& hs) const
 {
     dimensionedScalar Tstd("Tstd", dimTemperature, 298.15);
 
+    volScalarField htcst = htcs_->h()();
+    volScalarField htcwt = htcw_->h()();
+    forAll(alpha_, i)
+    {
+        htcst[i] *= max(alpha_[i], ROOTVSMALL);
+        htcwt[i] *= max(alpha_[i], ROOTVSMALL);
+    }
+    htcst.correctBoundaryConditions();
+    htcwt.correctBoundaryConditions();
+
     return
     (
-      - fvm::Sp(htcs_->h()/Cp_, hs) - htcs_->h()*(Tstd - TPrimary_)
-      - fvm::Sp(htcw_->h()/Cp_, hs) - htcw_->h()*(Tstd - Tw_)
+      - fvm::Sp(htcst/Cp_, hs) - htcst*(Tstd - TPrimary_)
+      - fvm::Sp(htcwt/Cp_, hs) -htcwt*(Tstd - Tw_)
     );
 }
 
diff --git a/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H b/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H
index c1269e8000..7ecbb075eb 100644
--- a/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H
+++ b/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,13 +35,17 @@ Description
 #include "specie.H"
 #include "rhoConst.H"
 #include "hConstThermo.H"
+#include "hPolynomialThermo.H"
 #include "hExponentialThermo.H"
 #include "constIsoSolidTransport.H"
 #include "constAnIsoSolidTransport.H"
 #include "exponentialSolidTransport.H"
+#include "polynomialSolidTransport.H"
+
 #include "sensibleEnthalpy.H"
 #include "thermo.H"
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -57,8 +61,7 @@ namespace Foam
                 >,
                 sensibleEnthalpy
             >
-        >
-        hConstSolidThermoPhysics;
+        > hConstSolidThermoPhysics;
 
     typedef
         exponentialSolidTransport
@@ -71,8 +74,23 @@ namespace Foam
                 >,
                 sensibleEnthalpy
             >
-        >
-        hExponentialSolidThermoPhysics;
+        > hExponentialSolidThermoPhysics;
+
+
+      typedef
+        polynomialSolidTransport
+        <
+            species::thermo
+            <
+                hPolynomialThermo
+                <
+                    rhoConst<specie>,
+                    8
+                >,
+                sensibleEnthalpy
+            >,
+            8
+        > hTransportThermoPoly8SolidThermoPhysics;
 }
 
 
diff --git a/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransport.C b/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransport.C
new file mode 100644
index 0000000000..7bddd6d822
--- /dev/null
+++ b/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransport.C
@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polynomialSolidTransport.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Thermo, int PolySize>
+Foam::polynomialSolidTransport<Thermo, PolySize>::polynomialSolidTransport
+(
+    Istream& is
+)
+:
+    Thermo(is),
+    kappaCoeffs_("kappaCoeffs<" + Foam::name(PolySize) + '>', is)
+{
+}
+
+
+template<class Thermo, int PolySize>
+Foam::polynomialSolidTransport<Thermo, PolySize>::polynomialSolidTransport
+(
+    const dictionary& dict
+)
+:
+    Thermo(dict),
+    kappaCoeffs_
+    (
+        dict.subDict("transport").lookup
+        (
+            "kappaCoeffs<" + Foam::name(PolySize) + '>'
+        )
+    )
+{
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Thermo, int PolySize>
+void Foam::polynomialSolidTransport<Thermo, PolySize>::write(Ostream& os) const
+{
+    os  << this->name() << endl;
+    os  << token::BEGIN_BLOCK << incrIndent << nl;
+
+    Thermo::write(os);
+
+    dictionary dict("transport");
+
+    dict.add
+    (
+        word("kappaCoeffs<" + Foam::name(PolySize) + '>'),
+        kappaCoeffs_
+    );
+    os  << indent << dict.dictName() << dict;
+
+    os  << decrIndent << token::END_BLOCK << nl;
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class Thermo, int PolySize>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const polynomialSolidTransport<Thermo, PolySize>& pt
+)
+{
+    os  << static_cast<const Thermo&>(pt) << tab
+        << "kappaCoeffs<" << Foam::name(PolySize) << '>' << tab
+        << pt.kappaCoeffs_;
+
+    os.check
+    (
+        "Ostream& operator<<"
+        "("
+            "Ostream&, "
+            "const polynomialSolidTransport<Thermo, PolySize>&"
+        ")"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransport.H b/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransport.H
new file mode 100644
index 0000000000..90a9e52f6a
--- /dev/null
+++ b/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransport.H
@@ -0,0 +1,237 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polynomialSolidTransport
+
+Description
+    Transport package using polynomial functions for solid kappa
+
+SourceFiles
+    polynomialSolidTransportI.H
+    polynomialSolidTransport.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polynomialSolidTransport_H
+#define polynomialSolidTransport_H
+
+#include "Polynomial.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Thermo, int PolySize> class polynomialSolidTransport;
+
+template<class Thermo, int PolySize>
+inline polynomialSolidTransport<Thermo, PolySize> operator+
+(
+    const polynomialSolidTransport<Thermo, PolySize>&,
+    const polynomialSolidTransport<Thermo, PolySize>&
+);
+
+template<class Thermo, int PolySize>
+inline polynomialSolidTransport<Thermo, PolySize> operator-
+(
+    const polynomialSolidTransport<Thermo, PolySize>&,
+    const polynomialSolidTransport<Thermo, PolySize>&
+);
+
+template<class Thermo, int PolySize>
+inline polynomialSolidTransport<Thermo, PolySize> operator*
+(
+    const scalar,
+    const polynomialSolidTransport<Thermo, PolySize>&
+);
+
+template<class Thermo, int PolySize>
+inline polynomialSolidTransport<Thermo, PolySize> operator==
+(
+    const polynomialSolidTransport<Thermo, PolySize>&,
+    const polynomialSolidTransport<Thermo, PolySize>&
+);
+
+template<class Thermo, int PolySize>
+Ostream& operator<<
+(
+    Ostream&,
+    const polynomialSolidTransport<Thermo, PolySize>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                     Class polynomialSolidTransport Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Thermo, int PolySize=8>
+class polynomialSolidTransport
+:
+    public Thermo
+{
+    // Private data
+
+        //- Thermal conductivity polynomial coefficients
+        //  Note: input in [W/m/K]
+        Polynomial<PolySize> kappaCoeffs_;
+
+
+    // Private Member Functions
+
+        //- Construct from components
+        inline polynomialSolidTransport
+        (
+            const Thermo& t,
+            const Polynomial<PolySize>& kappaPoly
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct copy
+        inline polynomialSolidTransport(const polynomialSolidTransport&);
+
+        //- Construct as named copy
+        inline polynomialSolidTransport
+        (
+            const word&,
+            const polynomialSolidTransport&
+        );
+
+        //- Construct from Istream
+        polynomialSolidTransport(Istream& is);
+
+        //- Construct from dictionary
+        polynomialSolidTransport(const dictionary& dict);
+
+        //- Construct and return a clone
+        inline autoPtr<polynomialSolidTransport> clone() const;
+
+        // Selector from Istream
+        inline static autoPtr<polynomialSolidTransport> New(Istream& is);
+
+        // Selector from dictionary
+        inline static autoPtr<polynomialSolidTransport> New
+        (
+            const dictionary&dict
+        );
+
+
+    // Member functions
+
+        //- Return the instantiated type name
+        static word typeName()
+        {
+            return "polynomial<" + Thermo::typeName() + '>';
+        }
+
+        //- Is the thermal conductivity isotropic
+        static const bool isotropic = true;
+
+        //- Dynamic viscosity [kg/ms]
+        inline scalar mu(const scalar p, const scalar T) const;
+
+        //- Thermal conductivity [W/mK]
+        inline scalar kappa(const scalar p, const scalar T) const;
+
+        //- Thermal conductivity [W/mK]
+        inline vector Kappa(const scalar p, const scalar T) const;
+
+        //- Thermal diffusivity of enthalpy [kg/ms]
+        inline scalar alphah(const scalar p, const scalar T) const;
+
+
+        //- Write to Ostream
+        void write(Ostream& os) const;
+
+
+    // Member operators
+
+        inline polynomialSolidTransport& operator=
+        (
+            const polynomialSolidTransport&
+        );
+        inline void operator+=(const polynomialSolidTransport&);
+        inline void operator-=(const polynomialSolidTransport&);
+        inline void operator*=(const scalar);
+
+
+    // Friend operators
+
+        friend polynomialSolidTransport operator+ <Thermo, PolySize>
+        (
+            const polynomialSolidTransport&,
+            const polynomialSolidTransport&
+        );
+
+        friend polynomialSolidTransport operator- <Thermo, PolySize>
+        (
+            const polynomialSolidTransport&,
+            const polynomialSolidTransport&
+        );
+
+        friend polynomialSolidTransport operator* <Thermo, PolySize>
+        (
+            const scalar,
+            const polynomialSolidTransport&
+        );
+
+        friend polynomialSolidTransport operator== <Thermo, PolySize>
+        (
+            const polynomialSolidTransport&,
+            const polynomialSolidTransport&
+        );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <Thermo, PolySize>
+        (
+            Ostream&,
+            const polynomialSolidTransport&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "polynomialSolidTransportI.H"
+
+#ifdef NoRepository
+#   include "polynomialSolidTransport.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransportI.H b/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransportI.H
new file mode 100644
index 0000000000..51dda6ba57
--- /dev/null
+++ b/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransportI.H
@@ -0,0 +1,289 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "specie.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Thermo, int PolySize>
+inline Foam::polynomialSolidTransport<Thermo, PolySize>::
+polynomialSolidTransport
+(
+    const polynomialSolidTransport& pt
+)
+:
+    Thermo(pt),
+    kappaCoeffs_(pt.kappaCoeffs_)
+{}
+
+
+template<class Thermo, int PolySize>
+inline Foam::polynomialSolidTransport<Thermo, PolySize>::
+polynomialSolidTransport
+(
+    const Thermo& t,
+    const Polynomial<PolySize>& kappaCoeffs
+)
+:
+    Thermo(t),
+    kappaCoeffs_(kappaCoeffs)
+{}
+
+
+template<class Thermo, int PolySize>
+inline Foam::polynomialSolidTransport<Thermo, PolySize>::
+polynomialSolidTransport
+(
+    const word& name,
+    const polynomialSolidTransport& pt
+)
+:
+    Thermo(name, pt),
+    kappaCoeffs_(pt.kappaCoeffs_)
+{}
+
+
+template<class Thermo, int PolySize>
+inline Foam::autoPtr<Foam::polynomialSolidTransport<Thermo, PolySize> >
+Foam::polynomialSolidTransport<Thermo, PolySize>::clone() const
+{
+    return autoPtr<polynomialSolidTransport<Thermo, PolySize> >
+    (
+        new polynomialSolidTransport<Thermo, PolySize>(*this)
+    );
+}
+
+
+template<class Thermo, int PolySize>
+inline Foam::autoPtr<Foam::polynomialSolidTransport<Thermo, PolySize> >
+Foam::polynomialSolidTransport<Thermo, PolySize>::New(Istream& is)
+{
+    return autoPtr<polynomialSolidTransport<Thermo, PolySize> >
+    (
+        new polynomialSolidTransport<Thermo, PolySize>(is)
+    );
+}
+
+
+template<class Thermo, int PolySize>
+inline Foam::autoPtr<Foam::polynomialSolidTransport<Thermo, PolySize> >
+Foam::polynomialSolidTransport<Thermo, PolySize>::New(const dictionary& dict)
+{
+    return autoPtr<polynomialSolidTransport<Thermo, PolySize> >
+    (
+        new polynomialSolidTransport<Thermo, PolySize>(dict)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Thermo, int PolySize>
+inline Foam::scalar Foam::polynomialSolidTransport<Thermo, PolySize>::mu
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    notImplemented
+    (
+        "Foam::scalar Foam::polynomialSolidTransport<thermo, PolySize>mu::"
+        "("
+        "    const scalar p, const scalar T"
+        ") const"
+    );
+    return scalar(0);
+}
+
+
+template<class Thermo, int PolySize>
+inline Foam::scalar Foam::polynomialSolidTransport<Thermo, PolySize>::kappa
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return kappaCoeffs_.value(T);
+}
+
+
+template<class Thermo, int PolySize>
+inline Foam::vector Foam::polynomialSolidTransport<Thermo, PolySize>::Kappa
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const scalar kappa(kappaCoeffs_.value(T));
+    return vector(kappa, kappa, kappa);
+}
+
+
+template<class Thermo, int PolySize>
+inline Foam::scalar Foam::polynomialSolidTransport<Thermo, PolySize>::alphah
+(
+    const scalar p, const scalar T
+) const
+{
+    return kappa(p, T)/this->Cpv(p, T);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Thermo, int PolySize>
+inline Foam::polynomialSolidTransport<Thermo, PolySize>&
+Foam::polynomialSolidTransport<Thermo, PolySize>::operator=
+(
+    const polynomialSolidTransport<Thermo, PolySize>& pt
+)
+{
+    Thermo::operator=(pt);
+
+    kappaCoeffs_ = pt.kappaCoeffs_;
+
+    return *this;
+}
+
+
+template<class Thermo, int PolySize>
+inline void Foam::polynomialSolidTransport<Thermo, PolySize>::operator+=
+(
+    const polynomialSolidTransport<Thermo, PolySize>& pt
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator+=(pt);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = pt.nMoles()/this->nMoles();
+
+    kappaCoeffs_ = molr1*kappaCoeffs_ + molr2*pt.kappaCoeffs_;
+}
+
+
+template<class Thermo, int PolySize>
+inline void Foam::polynomialSolidTransport<Thermo, PolySize>::operator-=
+(
+    const polynomialSolidTransport<Thermo, PolySize>& pt
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator-=(pt);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = pt.nMoles()/this->nMoles();
+
+    kappaCoeffs_ = molr1*kappaCoeffs_ - molr2*pt.kappaCoeffs_;
+}
+
+
+template<class Thermo, int PolySize>
+inline void Foam::polynomialSolidTransport<Thermo, PolySize>::operator*=
+(
+    const scalar s
+)
+{
+    Thermo::operator*=(s);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+template<class Thermo, int PolySize>
+inline Foam::polynomialSolidTransport<Thermo, PolySize> Foam::operator+
+(
+    const polynomialSolidTransport<Thermo, PolySize>& pt1,
+    const polynomialSolidTransport<Thermo, PolySize>& pt2
+)
+{
+    Thermo t
+    (
+        static_cast<const Thermo&>(pt1) + static_cast<const Thermo&>(pt2)
+    );
+
+    scalar molr1 = pt1.nMoles()/t.nMoles();
+    scalar molr2 = pt2.nMoles()/t.nMoles();
+
+    return polynomialSolidTransport<Thermo, PolySize>
+    (
+        t,
+        molr1*pt1.kappaCoeffs_ + molr2*pt2.kappaCoeffs_
+    );
+}
+
+
+template<class Thermo, int PolySize>
+inline Foam::polynomialSolidTransport<Thermo, PolySize> Foam::operator-
+(
+    const polynomialSolidTransport<Thermo, PolySize>& pt1,
+    const polynomialSolidTransport<Thermo, PolySize>& pt2
+)
+{
+    Thermo t
+    (
+        static_cast<const Thermo&>(pt1) - static_cast<const Thermo&>(pt2)
+    );
+
+    scalar molr1 = pt1.nMoles()/t.nMoles();
+    scalar molr2 = pt2.nMoles()/t.nMoles();
+
+    return polynomialSolidTransport<Thermo, PolySize>
+    (
+        t,
+        molr1*pt1.kappaCoeffs_ - molr2*pt2.kappaCoeffs_
+    );
+}
+
+
+template<class Thermo, int PolySize>
+inline Foam::polynomialSolidTransport<Thermo, PolySize> Foam::operator*
+(
+    const scalar s,
+    const polynomialSolidTransport<Thermo, PolySize>& pt
+)
+{
+    return polynomialSolidTransport<Thermo, PolySize>
+    (
+        s*static_cast<const Thermo&>(pt),
+        pt.kappaCoeffs_
+    );
+}
+
+
+template<class Thermo, int PolySize>
+inline Foam::polynomialSolidTransport<Thermo, PolySize> Foam::operator==
+(
+    const polynomialSolidTransport<Thermo, PolySize>& pt1,
+    const polynomialSolidTransport<Thermo, PolySize>& pt2
+)
+{
+    return pt2 - pt1;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
index 85000a3c10..5193a4b9de 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
+++ b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
@@ -103,6 +103,58 @@ addToRunTimeSelectionTable                                                    \
 );
 
 
+
+#define makeSolidThermoPhysicsType(BaseThermo,Cthermo,Mixture,SolidPhysicsType)\
+                                                                              \
+                                                                              \
+                                                                              \
+typedef                                                                       \
+    heThermo                                                                  \
+    <                                                                         \
+        BaseThermo,                                                           \
+        Mixture<SolidPhysicsType>                                             \
+    > heThermo##Mixture##SolidPhysicsType;                                    \
+                                                                              \
+typedef                                                                       \
+    Cthermo                                                                   \
+    <                                                                         \
+        BaseThermo,                                                           \
+        Mixture<SolidPhysicsType>                                             \
+    > Cthermo##Mixture##SolidPhysicsType;                                     \
+                                                                              \
+                                                                              \
+defineTemplateTypeNameAndDebugWithName                                        \
+(                                                                             \
+    Cthermo##Mixture##SolidPhysicsType,                                       \
+    (                                                                         \
+        #Cthermo"<"#Mixture"<"                                                \
+      + SolidPhysicsType::typeName()                                          \
+      + ">>"                                                                  \
+    ).c_str(),                                                                \
+    0                                                                         \
+);                                                                            \
+                                                                              \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicThermo,                                                              \
+    Cthermo##Mixture##SolidPhysicsType,                                       \
+    fvMesh                                                                    \
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    BaseThermo,                                                               \
+    Cthermo##Mixture##SolidPhysicsType,                                       \
+    fvMesh                                                                    \
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    BaseThermo,                                                               \
+    Cthermo##Mixture##SolidPhysicsType,                                       \
+    dictionary                                                                \
+);
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/solidThermos.C b/src/thermophysicalModels/solidThermo/solidThermo/solidThermos.C
index 6814431fec..020dbfda92 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/solidThermos.C
+++ b/src/thermophysicalModels/solidThermo/solidThermo/solidThermos.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,14 +31,17 @@ License
 #include "rhoConst.H"
 #include "hConstThermo.H"
 #include "hExponentialThermo.H"
+#include "hPolynomialThermo.H"
 #include "constIsoSolidTransport.H"
 #include "constAnIsoSolidTransport.H"
 #include "exponentialSolidTransport.H"
+#include "polynomialSolidTransport.H"
 #include "pureMixture.H"
 #include "sensibleEnthalpy.H"
 #include "sensibleInternalEnergy.H"
 #include "thermo.H"
 
+#include "solidThermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -83,6 +86,14 @@ makeSolidThermo
     specie
 );
 
+makeSolidThermoPhysicsType
+(
+    solidThermo,
+    heSolidThermo,
+    pureMixture,
+    hTransportThermoPoly8SolidThermoPhysics
+);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/thermophysicalModels/specie/transport/const/constTransportI.H b/src/thermophysicalModels/specie/transport/const/constTransportI.H
index 98ce21d97e..9949fddae4 100644
--- a/src/thermophysicalModels/specie/transport/const/constTransportI.H
+++ b/src/thermophysicalModels/specie/transport/const/constTransportI.H
@@ -111,7 +111,7 @@ inline Foam::scalar Foam::constTransport<Thermo>::kappa
     const scalar T
 ) const
 {
-    return this->Cp(p, T)*mu(p, T)*rPr_;
+    return this->Cpv(p, T)*mu(p, T)*rPr_;
 }
 
 
diff --git a/src/thermophysicalModels/specie/transport/polynomial/polynomialTransportI.H b/src/thermophysicalModels/specie/transport/polynomial/polynomialTransportI.H
index a2ff3f4d22..9a5ad80d5e 100644
--- a/src/thermophysicalModels/specie/transport/polynomial/polynomialTransportI.H
+++ b/src/thermophysicalModels/specie/transport/polynomial/polynomialTransportI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -129,7 +129,7 @@ inline Foam::scalar Foam::polynomialTransport<Thermo, PolySize>::alphah
     const scalar p, const scalar T
 ) const
 {
-    return kappa(p, T)/this->Cp(p, T);
+    return kappa(p, T)/this->Cpv(p, T);
 }
 
 
diff --git a/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H b/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H
index 7a0845cfb5..01668dda17 100644
--- a/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H
+++ b/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H
@@ -158,7 +158,7 @@ inline Foam::scalar Foam::sutherlandTransport<Thermo>::alphah
 ) const
 {
 
-    return kappa(p, T)/this->Cp(p, T);
+    return kappa(p, T)/this->Cpv(p, T);
 }
 
 
diff --git a/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.H b/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.H
index ba7fd1be49..e80d094844 100644
--- a/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.H
+++ b/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.H
@@ -80,7 +80,7 @@ public:
             //- Construct null
             inline deltaData();
 
-            //- Construct from origin, yStar, distance
+            //- Construct from delta value
             inline deltaData(const scalar delta);
 
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
index 5e78e48c97..3f1ab4a2c2 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
@@ -237,9 +237,9 @@ void epsilonWallFunctionFvPatchScalarField::calculate
 
         scalar w = cornerWeights[faceI];
 
-        epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
+        epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
 
-        G[cellI] =
+        G[cellI] +=
             w
            *(mutw[faceI] + muw[faceI])
            *magGradUw[faceI]
@@ -464,6 +464,8 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs
 
     FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());
 
+    scalarField& epsilonf = *this;
+
     // only set the values if the weights are < 1 - tolerance
     forAll(weights, faceI)
     {
@@ -475,6 +477,7 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs
 
             G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
             epsilon[cellI] = w*epsilon[cellI] + (1.0 - w)*epsilon0[cellI];
+            epsilonf[faceI] = epsilon[cellI];
         }
     }
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
index fff165e15e..2f4e11dc20 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -244,9 +244,9 @@ void omegaWallFunctionFvPatchScalarField::calculate
 
         scalar omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI]);
 
-        omega[cellI] = w*sqrt(sqr(omegaVis) + sqr(omegaLog));
+        omega[cellI] += w*sqrt(sqr(omegaVis) + sqr(omegaLog));
 
-        G[cellI] =
+        G[cellI] +=
             w
            *(mutw[faceI] + muw[faceI])
            *magGradUw[faceI]
@@ -481,6 +481,8 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs
 
     FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
 
+    scalarField& omegaf = *this;
+
     // only set the values if the weights are < 1 - tolerance
     forAll(weights, faceI)
     {
@@ -492,6 +494,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs
 
             G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
             omega[cellI] = w*omega[cellI] + (1.0 - w)*omega0[cellI];
+            omegaf[faceI] = omega[cellI];
         }
     }
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
index 7fb6fcfb96..dc08123d26 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
@@ -464,6 +464,8 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs
 
     FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());
 
+    scalarField& epsilonf = *this;
+
     // only set the values if the weights are < 1 - tolerance
     forAll(weights, faceI)
     {
@@ -475,6 +477,7 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs
 
             G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
             epsilon[cellI] = w*epsilon[cellI] + (1.0 - w)*epsilon0[cellI];
+            epsilonf[faceI] = epsilon[cellI];
         }
     }
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
index 554c18542d..69f8ffd170 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -242,9 +242,9 @@ void omegaWallFunctionFvPatchScalarField::calculate
 
         scalar omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI]);
 
-        omega[cellI] = w*sqrt(sqr(omegaVis) + sqr(omegaLog));
+        omega[cellI] += w*sqrt(sqr(omegaVis) + sqr(omegaLog));
 
-        G[cellI] =
+        G[cellI] +=
             w
            *(nutw[faceI] + nuw[faceI])
            *magGradUw[faceI]
@@ -479,6 +479,8 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs
 
     FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
 
+    scalarField& omegaf = *this;
+
     // only set the values if the weights are < 1 - tolerance
     forAll(weights, faceI)
     {
@@ -490,6 +492,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs
 
             G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
             omega[cellI] = w*omega[cellI] + (1.0 - w)*omega0[cellI];
+            omegaf[faceI] = omega[cellI];
         }
     }
 
diff --git a/wmake/wmake b/wmake/wmake
index ebcd38cb6e..b9977fa7e1 100755
--- a/wmake/wmake
+++ b/wmake/wmake
@@ -50,6 +50,7 @@ The 'target' is a Makefile target:
   e.g., Make/linux64GccDPOpt/fvMesh.o
 
 or a special target:
+  dwim              search up directories tree for Make sub-directory and build
   all               all subdirectories, uses any Allwmake files if they exist
   exe               build statically linked executable
   lib               build statically linked archive lib (.a)
@@ -136,7 +137,7 @@ fi
 
 
 #------------------------------------------------------------------------------
-# check arguments and change to the directory in which to run wmake
+# Check arguments and change to the directory in which to run wmake
 #------------------------------------------------------------------------------
 
 unset dir makeType
@@ -167,6 +168,51 @@ then
 fi
 
 
+#------------------------------------------------------------------------------
+# Search up directories tree for the Make sub-directory and build there
+#------------------------------------------------------------------------------
+
+unset dir
+MakeDir=Make
+
+expandpath()
+{
+    dir=`dirname $1`
+    cwd=$PWD
+    cd $dir
+    exPath=$PWD
+    cd $cwd
+}
+
+find_target()
+{
+    expandpath $1
+
+    if [ "$exPath" = "$WM_PROJECT_DIR" -o "$exPath" = "$HOME" -o "$exPath" = "/" ]; then
+        echo "$Script error: could not find Make directory"
+    elif [ -d "$1/Make" ]; then
+        echo Target $1
+        dir=$1
+    else
+        find_target "$1/.."
+    fi
+}
+
+if [ "$makeType" = dwim ]
+then
+    find_target .
+    makeType=
+
+    if [ "$dir" ]
+    then
+        cd $dir 2>/dev/null || {
+            echo "$Script error: could not change to directory '$dir'" 1>&2
+            exit 1
+        }
+    fi
+fi
+
+
 #------------------------------------------------------------------------------
 # Recurse the application directories tree
 #------------------------------------------------------------------------------
@@ -254,8 +300,9 @@ $make -s -f $WM_DIR/Makefile MAKE_DIR=$MakeDir INCLUDE_DEPS=$OBJECTS_DIR/dontInc
 rc=$?
 [ $rc -eq 0 ] || exit $rc
 
+
 #------------------------------------------------------------------------------
-# make the object files and link
+# Make the object files and link
 #------------------------------------------------------------------------------
 
 cmd="$make -f $WM_DIR/Makefile MAKE_DIR=$MakeDir INCLUDE_DEPS=$OBJECTS_DIR/includeDeps $makeType"