diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H
index 1fe150e1..2a5f6c8b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H
@@ -14,7 +14,6 @@
 
       + fvc::ddt(alpha1, rho1, K1) + fvc::div(alphaRhoPhi1, K1)
       - contErr1*K1
-
       + (
             he1.name() == thermo1.phasePropertyName("e")
           ? fvc::ddt(alpha1)*p + fvc::div(alphaPhi1, p)
@@ -27,9 +26,12 @@
            *fvc::interpolate(thermo1.alphaEff(phase1.turbulence().mut())),
             he1
         )
+    );
 
-     ==
+    he1Eqn.relax();
 
+    he1Eqn -=
+    (
         heatTransferCoeff*(thermo2.T() - thermo1.T())
       + heatTransferCoeff*he1/Cpv1
       - fvm::Sp(heatTransferCoeff/Cpv1, he1)
@@ -43,7 +45,6 @@
 
       + fvc::ddt(alpha2, rho2, K2) + fvc::div(alphaRhoPhi2, K2)
       - contErr2*K2
-
       + (
             he2.name() == thermo2.phasePropertyName("e")
           ? fvc::ddt(alpha2)*p + fvc::div(alphaPhi2, p)
@@ -56,23 +57,29 @@
            *fvc::interpolate(thermo2.alphaEff(phase2.turbulence().mut())),
             he2
         )
+    );
 
-     ==
+    he2Eqn.relax();
 
+    he2Eqn -=
+    (
         heatTransferCoeff*(thermo1.T() - thermo2.T())
       + heatTransferCoeff*he2/Cpv2
       - fvm::Sp(heatTransferCoeff/Cpv2, he2)
       + fvOptions(alpha2, rho2, he2)
     );
 
-    he1Eqn.relax();
     fvOptions.constrain(he1Eqn);
     he1Eqn.solve();
 
-    he2Eqn.relax();
     fvOptions.constrain(he2Eqn);
     he2Eqn.solve();
 
     thermo1.correct();
+    Info<< "min " << thermo1.T().name()
+        << " " << min(thermo1.T()).value() << endl;
+
     thermo2.correct();
+    Info<< "min " << thermo2.T().name()
+        << " " << min(thermo2.T()).value() << endl;
 }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
index 45a52bc5..bfd5faf8 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
@@ -16,13 +16,11 @@ volScalarField dragCoeff(fluid.dragCoeff());
     {
         U1Eqn =
         (
-            fvm::ddt(alpha1, rho1, U1)
-          + fvm::div(alphaRhoPhi1, U1)
+            fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
           - fvm::Sp(contErr1, U1)
           + mrfZones(alpha1*rho1 + virtualMassCoeff, U1)
           + phase1.turbulence().divDevRhoReff(U1)
          ==
-          - fvm::Sp(dragCoeff, U1)
           - liftForce
           - wallLubricationForce
           - turbulentDispersionForce
@@ -36,20 +34,19 @@ volScalarField dragCoeff(fluid.dragCoeff());
           + fvOptions(alpha1, rho1, U1)
         );
         U1Eqn.relax();
+        U1Eqn += fvm::Sp(dragCoeff, U1);
         fvOptions.constrain(U1Eqn);
     }
 
     {
         U2Eqn =
         (
-            fvm::ddt(alpha2, rho2, U2)
-          + fvm::div(alphaRhoPhi2, U2)
+            fvm::ddt(alpha2, rho2, U2) + fvm::div(alphaRhoPhi2, U2)
           - fvm::Sp(contErr2, U2)
           + mrfZones(alpha2*rho2 + virtualMassCoeff, U2)
           + phase2.turbulence().divDevRhoReff(U2)
          ==
-          - fvm::Sp(dragCoeff, U2)
-          + liftForce
+            liftForce
           + wallLubricationForce
           + turbulentDispersionForce
           - virtualMassCoeff
@@ -62,6 +59,7 @@ volScalarField dragCoeff(fluid.dragCoeff());
           + fvOptions(alpha2, rho2, U2)
         );
         U2Eqn.relax();
+        U2Eqn += fvm::Sp(dragCoeff, U2);
         fvOptions.constrain(U2Eqn);
     }
 }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/files b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/files
index 81bccf25..ce1bdbc4 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/files
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/files
@@ -56,4 +56,6 @@ aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.C
 aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C
 aspectRatioModels/Wellek/Wellek.C
 
+wallDependentModel/wallDependentModel.C
+
 LIB = $(FOAM_LIBBIN)/libcompressibleEulerianInterfacialModels
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/options b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/options
index a728cb7f..25d20582 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/options
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/transportModel \
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.C
index cd254865..7660bb97 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.C
@@ -53,7 +53,7 @@ Foam::aspectRatioModels::TomiyamaAspectRatio::TomiyamaAspectRatio
 )
 :
     VakhrushevEfremov(dict, pair),
-    yWall_(pair.phase1().mesh().lookupObject<volScalarField>("yWall"))
+    wallDependentModel(pair.phase1().mesh())
 {}
 
 
@@ -72,7 +72,7 @@ Foam::aspectRatioModels::TomiyamaAspectRatio::E() const
         VakhrushevEfremov::E()
        *max
        (
-           scalar(1) - 0.35*yWall_/pair_.dispersed().d(),
+           scalar(1) - 0.35*yWall()/pair_.dispersed().d(),
            scalar(0.65)
        );
 }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.H
index 82f00642..747214f9 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.H
@@ -45,6 +45,7 @@ SourceFiles
 #define TomiyamaAspectRatio_H
 
 #include "VakhrushevEfremov.H"
+#include "wallDependentModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,16 +60,9 @@ namespace aspectRatioModels
 
 class TomiyamaAspectRatio
 :
-    public VakhrushevEfremov
+    public VakhrushevEfremov,
+    public wallDependentModel
 {
-private:
-
-    // Private data
-
-        //- Wall distance
-        const volScalarField& yWall_;
-
-
 public:
 
     //- Runtime type information
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallDependentModel/wallDependentModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallDependentModel/wallDependentModel.C
new file mode 100644
index 00000000..b64547a1
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallDependentModel/wallDependentModel.C
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenCFD Ltd.
+     \\/     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 "wallDependentModel.H"
+#include "wallDist.H"
+#include "wallDistReflection.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallDependentModel::wallDependentModel(const fvMesh& mesh)
+:
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::wallDependentModel::~wallDependentModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const Foam::volScalarField& Foam::wallDependentModel::yWall() const
+{
+    if (!mesh_.foundObject<volScalarField>("yWall"))
+    {
+        wallDist w(mesh_);
+
+        volScalarField* yPtr
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "yWall",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    true
+                ),
+                w.y()
+            )
+        );
+
+        yPtr->checkIn();
+    }
+
+    return mesh_.lookupObject<volScalarField>("yWall");
+}
+
+
+const Foam::volVectorField& Foam::wallDependentModel::nWall() const
+{
+    if (!mesh_.foundObject<volVectorField>("nWall"))
+    {
+        wallDistReflection w(mesh_);
+
+        if (!mesh_.foundObject<volScalarField>("yWall"))
+        {
+            volScalarField* yPtr
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "yWall",
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE,
+                        true
+                    ),
+                    w.y()
+                )
+            );
+
+            yPtr->checkIn();
+        }
+
+        volVectorField* nPtr
+        (
+            new volVectorField
+            (
+                IOobject
+                (
+                    "nWall",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    true
+                ),
+                w.n()
+            )
+        );
+
+        nPtr->checkIn();
+    }
+
+    return mesh_.lookupObject<volVectorField>("nWall");
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallDependentModel/wallDependentModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallDependentModel/wallDependentModel.H
new file mode 100644
index 00000000..411add08
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallDependentModel/wallDependentModel.H
@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenCFD Ltd.
+     \\/     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::wallDependentModel
+
+Description
+    A class which provides on-demand creation and caching of wall distance and
+    wall normal fields for use by multiple models.
+
+SourceFiles
+    wallDependentModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallDependentModel_H
+#define wallDependentModel_H
+
+#include "fvMesh.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class wallDependentModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallDependentModel
+{
+    // Private data
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        wallDependentModel(const wallDependentModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const wallDependentModel&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from a mesh
+        wallDependentModel(const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~wallDependentModel();
+
+    // Member Functions
+
+        // Return the wall distance, creating and storing it if necessary
+        const volScalarField& yWall() const;
+
+        // Return the wall normal, creating and storing it if necessary
+        const volVectorField& nWall() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.C
index 9b047289..b2568dd7 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.C
@@ -70,19 +70,19 @@ Foam::wallLubricationModels::Antal::~Antal()
 Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Antal::F() const
 {
     volVectorField Ur(pair_.Ur());
-    volVectorField nWall(- fvc::grad(yWall_));
-    nWall /= mag(nWall) + SMALL;
+
+    const volVectorField& n(nWall());
 
     return
         max
         (
             dimensionedScalar("zero", dimless/dimLength, 0),
-            Cw1_/pair_.dispersed().d() + Cw2_/yWall_
+            Cw1_/pair_.dispersed().d() + Cw2_/yWall()
         )
        *pair_.dispersed()
        *pair_.continuous().rho()
-       *magSqr(Ur - (Ur & nWall)*nWall)
-       *nWall;
+       *magSqr(Ur - (Ur & n)*n)
+       *n;
 }
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C
index c89f1233..1ec5df33 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C
@@ -71,11 +71,12 @@ Foam::wallLubricationModels::Frank::~Frank()
 Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Frank::F() const
 {
     volVectorField Ur(pair_.Ur());
-    volVectorField nWall(- fvc::grad(yWall_));
-    nWall /= mag(nWall) + SMALL;
+
+    const volVectorField& n(nWall());
+    const volScalarField& y(yWall());
 
     volScalarField Eo(pair_.Eo());
-    volScalarField yTilde(yWall_/(Cwc_*pair_.dispersed().d()));
+    volScalarField yTilde(y/(Cwc_*pair_.dispersed().d()));
 
     return
         (
@@ -86,12 +87,12 @@ Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Frank::F() const
        *max
         (
             dimensionedScalar("zero", dimless/dimLength, 0.0),
-            (1.0 - yTilde)/(Cwd_*yWall_*pow(yTilde, p_ - 1.0))
+            (1.0 - yTilde)/(Cwd_*y*pow(yTilde, p_ - 1.0))
         )
        *pair_.dispersed()
        *pair_.continuous().rho()
-       *magSqr(Ur - (Ur & nWall)*nWall)
-       *nWall;
+       *magSqr(Ur - (Ur & n)*n)
+       *n;
 }
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C
index 09e45c32..b220aab6 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C
@@ -70,8 +70,9 @@ Foam::tmp<Foam::volVectorField>
 Foam::wallLubricationModels::TomiyamaWallLubrication::F() const
 {
     volVectorField Ur(pair_.Ur());
-    volVectorField nWall(- fvc::grad(yWall_));
-    nWall /= mag(nWall) + SMALL;
+
+    const volVectorField& n(nWall());
+    const volScalarField& y(yWall());
 
     volScalarField Eo(pair_.Eo());
 
@@ -84,13 +85,13 @@ Foam::wallLubricationModels::TomiyamaWallLubrication::F() const
        *0.5
        *pair_.dispersed().d()
        *(
-            1/sqr(yWall_)
-          - 1/sqr(D_ - yWall_)
+            1/sqr(y)
+          - 1/sqr(D_ - y)
         )
        *pair_.dispersed()
        *pair_.continuous().rho()
-       *magSqr(Ur - (Ur & nWall)*nWall)
-       *nWall;
+       *magSqr(Ur - (Ur & n)*n)
+       *n;
 }
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
index 593516e4..b7a7edc3 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
@@ -45,8 +45,8 @@ Foam::wallLubricationModel::wallLubricationModel
     const phasePair& pair
 )
 :
-    pair_(pair),
-    yWall_(pair.phase1().mesh().lookupObject<volScalarField>("yWall"))
+    wallDependentModel(pair.phase1().mesh()),
+    pair_(pair)
 {}
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.H
index 04fd4729..8f77e213 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.H
@@ -37,6 +37,7 @@ SourceFiles
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "wallDependentModel.H"
 #include "volFields.H"
 #include "dictionary.H"
 #include "runTimeSelectionTables.H"
@@ -51,6 +52,8 @@ class phasePair;
 \*---------------------------------------------------------------------------*/
 
 class wallLubricationModel
+:
+    public wallDependentModel
 {
 protected:
 
@@ -59,9 +62,6 @@ protected:
         //- Phase pair
         const phasePair& pair_;
 
-        //- Wall distance
-        const volScalarField& yWall_;
-
 
 public:
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index f8a18f4f..dcc3595b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -32,7 +32,6 @@ License
 #include "liftModel.H"
 #include "wallLubricationModel.H"
 #include "turbulentDispersionModel.H"
-#include "wallDist.H"
 #include "fvMatrix.H"
 #include "surfaceInterpolate.H"
 #include "MULES.H"
@@ -110,17 +109,6 @@ Foam::twoPhaseSystem::twoPhaseSystem
         ),
         mesh,
         dimensionedScalar("dgdt", dimless/dimTime, 0)
-    ),
-
-    yWall_
-    (
-        IOobject
-        (
-            "yWall",
-            mesh.time().timeName(),
-            mesh
-        ),
-        wallDist(mesh).y()
     )
 {
     phase2_.volScalarField::operator=(scalar(1) - phase1_);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
index e6a5732b..a12f3a4e 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
@@ -83,9 +83,6 @@ private:
         //-  Dilatation term
         volScalarField dgdt_;
 
-        //- Wall distance
-        volScalarField yWall_;
-
         //- Unordered phase pair
         autoPtr<phasePair> pair_;
 
diff --git a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C
index 1f91954e..dc797fdc 100644
--- a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C
+++ b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C
@@ -53,7 +53,8 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
-#   include "addOverwriteOption.H"
+    #include "addOverwriteOption.H"
+    #include "addRegionOption.H"
     argList::validArgs.append("cellSet");
     argList::addBoolOption
     (
@@ -62,10 +63,10 @@ int main(int argc, char *argv[])
         " (default is to extend set)"
     );
 
-#   include "setRootCase.H"
-#   include "createTime.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
     runTime.functionObjects().off();
-#   include "createMesh.H"
+    #include "createNamedMesh.H"
     const word oldInstance = mesh.pointsInstance();
 
     word cellSetName(args.args()[1]);
@@ -181,7 +182,7 @@ int main(int argc, char *argv[])
         mesh.movePoints(map().preMotionPoints());
     }
 
-    Pout<< "Refined from " << returnReduce(map().nOldCells(), sumOp<label>())
+    Info<< "Refined from " << returnReduce(map().nOldCells(), sumOp<label>())
         << " to " << mesh.globalData().nTotalCells() << " cells." << nl << endl;
 
     if (overwrite)
diff --git a/applications/utilities/miscellaneous/patchSummary/patchSummary.C b/applications/utilities/miscellaneous/patchSummary/patchSummary.C
index e26fc109..01b2137f 100644
--- a/applications/utilities/miscellaneous/patchSummary/patchSummary.C
+++ b/applications/utilities/miscellaneous/patchSummary/patchSummary.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,7 @@ Description
 
 #include "fvCFD.H"
 #include "volFields.H"
+#include "pointFields.H"
 #include "IOobjectList.H"
 #include "patchSummaryTemplates.H"
 
@@ -87,6 +88,12 @@ int main(int argc, char *argv[])
         PtrList<volSymmTensorField> vsytf(objNames.size());
         PtrList<volTensorField> vtf(objNames.size());
 
+        PtrList<pointScalarField> psf(objNames.size());
+        PtrList<pointVectorField> pvf(objNames.size());
+        PtrList<pointSphericalTensorField> psptf(objNames.size());
+        PtrList<pointSymmTensorField> psytf(objNames.size());
+        PtrList<pointTensorField> ptf(objNames.size());
+
         Info<< "Valid fields:" << endl;
 
         forAll(objNames, objI)
@@ -101,11 +108,17 @@ int main(int argc, char *argv[])
 
             if (obj.headerOk())
             {
-                addToFieldList<scalar>(vsf, obj, objI, mesh);
-                addToFieldList<vector>(vvf, obj, objI, mesh);
-                addToFieldList<sphericalTensor>(vsptf, obj, objI, mesh);
-                addToFieldList<symmTensor>(vsytf, obj, objI, mesh);
-                addToFieldList<tensor>(vtf, obj, objI, mesh);
+                addToFieldList(vsf, obj, objI, mesh);
+                addToFieldList(vvf, obj, objI, mesh);
+                addToFieldList(vsptf, obj, objI, mesh);
+                addToFieldList(vsytf, obj, objI, mesh);
+                addToFieldList(vtf, obj, objI, mesh);
+
+                addToFieldList(psf, obj, objI, pointMesh::New(mesh));
+                addToFieldList(pvf, obj, objI, pointMesh::New(mesh));
+                addToFieldList(psptf, obj, objI, pointMesh::New(mesh));
+                addToFieldList(psytf, obj, objI, pointMesh::New(mesh));
+                addToFieldList(ptf, obj, objI, pointMesh::New(mesh));
             }
         }
 
@@ -119,11 +132,17 @@ int main(int argc, char *argv[])
             forAll(bm, patchI)
             {
                 Info<< bm[patchI].type() << "\t: " << bm[patchI].name() << nl;
-                outputFieldList<scalar>(vsf, patchI);
-                outputFieldList<vector>(vvf, patchI);
-                outputFieldList<sphericalTensor>(vsptf, patchI);
-                outputFieldList<symmTensor>(vsytf, patchI);
-                outputFieldList<tensor>(vtf, patchI);
+                outputFieldList(vsf, patchI);
+                outputFieldList(vvf, patchI);
+                outputFieldList(vsptf, patchI);
+                outputFieldList(vsytf, patchI);
+                outputFieldList(vtf, patchI);
+
+                outputFieldList(psf, patchI);
+                outputFieldList(pvf, patchI);
+                outputFieldList(psptf, patchI);
+                outputFieldList(psytf, patchI);
+                outputFieldList(ptf, patchI);
                 Info<< endl;
             }
         }
@@ -139,11 +158,17 @@ int main(int argc, char *argv[])
             forAll(bm, patchI)
             {
                 HashTable<word> fieldToType;
-                collectFieldList<scalar>(vsf, patchI, fieldToType);
-                collectFieldList<vector>(vvf, patchI, fieldToType);
-                collectFieldList<sphericalTensor>(vsptf, patchI, fieldToType);
-                collectFieldList<symmTensor>(vsytf, patchI, fieldToType);
-                collectFieldList<tensor>(vtf, patchI, fieldToType);
+                collectFieldList(vsf, patchI, fieldToType);
+                collectFieldList(vvf, patchI, fieldToType);
+                collectFieldList(vsptf, patchI, fieldToType);
+                collectFieldList(vsytf, patchI, fieldToType);
+                collectFieldList(vtf, patchI, fieldToType);
+
+                collectFieldList(psf, patchI, fieldToType);
+                collectFieldList(pvf, patchI, fieldToType);
+                collectFieldList(psptf, patchI, fieldToType);
+                collectFieldList(psytf, patchI, fieldToType);
+                collectFieldList(ptf, patchI, fieldToType);
 
                 label groupI = findIndex(fieldToTypes, fieldToType);
                 if (groupI == -1)
@@ -184,11 +209,17 @@ int main(int argc, char *argv[])
                             Info<< "group\t: " << groups[i] << nl;
                         }
                     }
-                    outputFieldList<scalar>(vsf, patchIDs[0]);
-                    outputFieldList<vector>(vvf, patchIDs[0]);
-                    outputFieldList<sphericalTensor>(vsptf, patchIDs[0]);
-                    outputFieldList<symmTensor>(vsytf, patchIDs[0]);
-                    outputFieldList<tensor>(vtf, patchIDs[0]);
+                    outputFieldList(vsf, patchIDs[0]);
+                    outputFieldList(vvf, patchIDs[0]);
+                    outputFieldList(vsptf, patchIDs[0]);
+                    outputFieldList(vsytf, patchIDs[0]);
+                    outputFieldList(vtf, patchIDs[0]);
+
+                    outputFieldList(psf, patchIDs[0]);
+                    outputFieldList(pvf, patchIDs[0]);
+                    outputFieldList(psptf, patchIDs[0]);
+                    outputFieldList(psytf, patchIDs[0]);
+                    outputFieldList(ptf, patchIDs[0]);
                     Info<< endl;
                 }
                 else
@@ -199,11 +230,17 @@ int main(int argc, char *argv[])
                         label patchI = patchIDs[i];
                         Info<< bm[patchI].type()
                             << "\t: " << bm[patchI].name() << nl;
-                        outputFieldList<scalar>(vsf, patchI);
-                        outputFieldList<vector>(vvf, patchI);
-                        outputFieldList<sphericalTensor>(vsptf, patchI);
-                        outputFieldList<symmTensor>(vsytf, patchI);
-                        outputFieldList<tensor>(vtf, patchI);
+                        outputFieldList(vsf, patchI);
+                        outputFieldList(vvf, patchI);
+                        outputFieldList(vsptf, patchI);
+                        outputFieldList(vsytf, patchI);
+                        outputFieldList(vtf, patchI);
+
+                        outputFieldList(psf, patchI);
+                        outputFieldList(pvf, patchI);
+                        outputFieldList(psptf, patchI);
+                        outputFieldList(psytf, patchI);
+                        outputFieldList(ptf, patchI);
                         Info<< endl;
                     }
                 }
diff --git a/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.C b/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.C
index bc0eee1e..4f130eed 100644
--- a/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.C
+++ b/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,33 +28,33 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class Type>
+template<class GeoField>
 void Foam::addToFieldList
 (
-    PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+    PtrList<GeoField>& fieldList,
     const IOobject& obj,
     const label fieldI,
-    const fvMesh& mesh
+    const typename GeoField::Mesh& mesh
 )
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+    typedef GeoField fieldType;
 
-    if (obj.headerClassName() == fieldType::typeName)
+    if (obj.headerClassName() == GeoField::typeName)
     {
         fieldList.set
         (
             fieldI,
-            new fieldType(obj, mesh)
+            new GeoField(obj, mesh)
         );
-        Info<< "    " << fieldType::typeName << tab << obj.name() << endl;
+        Info<< "    " << GeoField::typeName << tab << obj.name() << endl;
     }
 }
 
 
-template<class Type>
+template<class GeoField>
 void Foam::outputFieldList
 (
-    const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+    const PtrList<GeoField>& fieldList,
     const label patchI
 )
 {
@@ -62,7 +62,8 @@ void Foam::outputFieldList
     {
         if (fieldList.set(fieldI))
         {
-            Info<< "    " << pTraits<Type>::typeName << tab << tab
+            Info<< "    " << pTraits<typename GeoField::value_type>::typeName
+                << tab << tab
                 << fieldList[fieldI].name() << tab << tab
                 << fieldList[fieldI].boundaryField()[patchI].type() << nl;
         }
@@ -70,10 +71,10 @@ void Foam::outputFieldList
 }
 
 
-template<class Type>
+template<class GeoField>
 void Foam::collectFieldList
 (
-    const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+    const PtrList<GeoField>& fieldList,
     const label patchI,
     HashTable<word>& fieldToType
 )
diff --git a/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.H b/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.H
index d61a416a..8689ed84 100644
--- a/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.H
+++ b/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,26 +33,26 @@ License
 
 namespace Foam
 {
-    template<class Type>
+    template<class GeoField>
     void addToFieldList
     (
-        PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+        PtrList<GeoField>& fieldList,
         const IOobject& obj,
         const label fieldI,
-        const fvMesh& mesh
+        const typename GeoField::Mesh& mesh
     );
 
-    template<class Type>
+    template<class GeoField>
     void outputFieldList
     (
-        const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+        const PtrList<GeoField>& fieldList,
         const label patchI
     );
 
-    template<class Type>
+    template<class GeoField>
     void collectFieldList
     (
-        const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+        const PtrList<GeoField>& fieldList,
         const label patchI,
         HashTable<word>& fieldToType
     );
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/inletOutletTotalTemperature/inletOutletTotalTemperatureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/inletOutletTotalTemperature/inletOutletTotalTemperatureFvPatchScalarField.C
index a1c7e2ad..607f8893 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/inletOutletTotalTemperature/inletOutletTotalTemperatureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/inletOutletTotalTemperature/inletOutletTotalTemperatureFvPatchScalarField.C
@@ -77,7 +77,7 @@ inletOutletTotalTemperatureFvPatchScalarField
 :
     inletOutletFvPatchScalarField(p, iF),
     UName_(dict.lookupOrDefault<word>("U", "U")),
-    psiName_(dict.lookupOrDefault<word>("psi", "psi")),
+    psiName_(dict.lookupOrDefault<word>("psi", "thermo:psi")),
     gamma_(readScalar(dict.lookup("gamma"))),
     T0_("T0", dict, p.size())
 {
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
index 37c98e6a..871c7f8a 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -1177,6 +1177,15 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
         {
             newPoints[meshPoints[i]] += disp[i];
         }
+
+        syncTools::syncPointList
+        (
+            mesh_,
+            newPoints,
+            minMagSqrEqOp<point>(),     // combine op
+            vector(GREAT, GREAT, GREAT) // null value (note: cannot use VGREAT)
+        );
+
         mesh_.movePoints(newPoints);
     }
 
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C
index 59d8c5cd..4f4f8e5a 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C
@@ -97,11 +97,13 @@ void Foam::mapNearestAMI<SourcePatch, TargetPatch>::setNextNearestFaces
         }
     }
 
-    forAll(mapFlag, srcFaceI)
+    forAll(mapFlag, faceI)
     {
-        if (!mapFlag[srcFaceI])
+        if (mapFlag[faceI])
         {
-            tgtFaceI = this->findTargetFace(srcFaceI);
+            srcFaceI = faceI;
+            tgtFaceI = this->findTargetFace(faceI);
+
             if (tgtFaceI == -1)
             {
                 const vectorField& srcCf = this->srcPatch_.faceCentres();