diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
index 3fb6dfe45e412afeca48b498fb99a2127b03196a..7a5c717cff712da15a12e99911e4580fe6a96845 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
@@ -27,7 +27,7 @@
         fvScalarMatrix pDDtEqn
         (
             fvc::ddt(rho) + fvc::div(phiHbyA)
-          + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
+          + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
         );
 
         while (pimple.correctNonOrthogonal())
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
index 117790446e9509ee3dcbfd235f9dd8db44508a41..9bed803d1e7be78b71e1d7bfde39bb93f8bf1d0e 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
@@ -68,4 +68,6 @@
     }
 }
 
+phi.oldTime() = phi;
+
 #include "continuityErrs.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
index 92a2bd64d8f5d5b1317e54f80370ae06d87e46b1..d97e8b2a3549db8ea23edc6c073a5bb160c6537c 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
@@ -1,11 +1,11 @@
 {
-    solve
+    fvScalarMatrix TEqn
     (
         fvm::ddt(rho, T)
       + fvm::div(rhoPhi, T)
       - fvm::laplacian(twoPhaseProperties.alphaEff(turbulence->mut()), T)
       + (
-            p*fvc::div(phi)
+            fvc::div(fvc::absolute(phi, U), p)
           + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
         )
        *(
@@ -14,5 +14,8 @@
         )
     );
 
+    TEqn.relax();
+    TEqn.solve();
+
     twoPhaseProperties.correct();
 }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
index 1d7b9ca624a4e0acb1d653890c23adda345b5627..513ef961bd7e411cde9fee70434c1c418e84f4e5 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
@@ -56,6 +56,8 @@
             phi -= pcorrEqn.flux();
         }
     }
-
-    #include "continuityErrs.H"
 }
+
+phi.oldTime() = phi;
+
+#include "continuityErrs.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index 3148382ad2e23cb2ddc774dd7064b8bc28ce8c38..a42d9119804c8e61031ee217ec38b8722781f236 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -104,8 +104,6 @@ int main(int argc, char *argv[])
             }
         }
 
-        rho = alpha1*rho1 + alpha2*rho2;
-
         runTime.write();
 
         Info<< "ExecutionTime = "
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 102d4d57a99a8bfa111eaca227399c27d2066935..c7289b23f901906c4fe7c96366c041ac8f0a3e9d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -38,9 +38,9 @@
 
     volScalarField& p = twoPhaseProperties.p();
     volScalarField& T = twoPhaseProperties.T();
-    const volScalarField& rho1 = twoPhaseProperties.thermo1().rho();
+    volScalarField& rho1 = twoPhaseProperties.thermo1().rho();
     const volScalarField& psi1 = twoPhaseProperties.thermo1().psi();
-    const volScalarField& rho2 = twoPhaseProperties.thermo2().rho();
+    volScalarField& rho2 = twoPhaseProperties.thermo2().rho();
     const volScalarField& psi2 = twoPhaseProperties.thermo2().psi();
 
     volScalarField rho
@@ -93,18 +93,5 @@
         compressible::turbulenceModel::New(rho, U, rhoPhi, twoPhaseProperties)
     );
 
-    Info<< "Creating field dpdt\n" << endl;
-    volScalarField dpdt
-    (
-        IOobject
-        (
-            "dpdt",
-            runTime.timeName(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-    );
-
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 71745c0168b7bd70b0e8aa6d8ca376c92ab0bdfd..73babb08f04ae70d23f2b3ed58831060ec0c771e 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -26,28 +26,44 @@
     tmp<fvScalarMatrix> p_rghEqnComp1;
     tmp<fvScalarMatrix> p_rghEqnComp2;
 
-    //if (transonic)
-    //{
-    //}
-    //else
+    if (pimple.transonic())
     {
         surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
         surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
 
+        p_rghEqnComp1 =
+            fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
+          + correction
+            (
+                psi1*fvm::ddt(p_rgh)
+              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+            );
+        deleteDemandDrivenData(p_rghEqnComp1().faceFluxCorrectionPtr());
+        p_rghEqnComp1().relax();
+
+        p_rghEqnComp2 =
+            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+          + correction
+            (
+                psi2*fvm::ddt(p_rgh)
+              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+            );
+        deleteDemandDrivenData(p_rghEqnComp2().faceFluxCorrectionPtr());
+        p_rghEqnComp2().relax();
+    }
+    else
+    {
         p_rghEqnComp1 =
             fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
-          + fvc::div(phid1, p_rgh)
-          - fvc::Sp(fvc::div(phid1), p_rgh);
+          + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
 
         p_rghEqnComp2 =
             fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
-          + fvc::div(phid2, p_rgh)
-          - fvc::Sp(fvc::div(phid2), p_rgh);
+          + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
     }
 
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    //twoPhaseProperties.rho() -= psi*p_rgh;
+    // Cache p_rgh prior to solve for density update
+    volScalarField p_rgh_0(p_rgh);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -69,8 +85,8 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            // Second part of thermodynamic density update
-            //twoPhaseProperties.rho() += psi*p_rgh;
+            //p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+            //p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
 
             dgdt =
             (
@@ -88,15 +104,14 @@
 
     p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
 
-    // twoPhaseProperties.correct();
+    // Update densities from change in p_rgh
+    rho1 += psi1*(p_rgh - p_rgh_0);
+    rho2 += psi2*(p_rgh - p_rgh_0);
 
-    Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
+    rho = alpha1*rho1 + alpha2*rho2;
 
     K = 0.5*magSqr(U);
 
-    if (twoPhaseProperties.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
+    Info<< "max(U) " << max(mag(U)).value() << endl;
+    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
 }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H
index b2e636a4c471ac92b490aa58f1084b72c6dc3888..608dc91a9b18718a506dbb403a83203a4f6c9d87 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H
@@ -94,6 +94,16 @@ public:
             return thermo2_();
         }
 
+        rhoThermo& thermo1()
+        {
+            return thermo1_();
+        }
+
+        rhoThermo& thermo2()
+        {
+            return thermo2_();
+        }
+
         //- Update properties
         virtual void correct();
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
index 54a61f51d2c386470da4a91a9f687892ec8dc334..cc5b3ebe21915bc2b19c97c32a3bc067a28b73d7 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
@@ -62,10 +62,10 @@
 
     volScalarField& p = thermo1.p();
 
-    volScalarField rho1("rho" + phase1Name, thermo1.rho());
+    volScalarField& rho1 = thermo1.rho();
     const volScalarField& psi1 = thermo1.psi();
 
-    volScalarField rho2("rho" + phase2Name, thermo2.rho());
+    volScalarField& rho2 = thermo2.rho();
     const volScalarField& psi2 = thermo2.psi();
 
     volVectorField U
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index bf75b28d058bedc8da84271efc9a08de92d18f7d..02e49b56618430559a3482700ec6fa02971f22c9 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -1,7 +1,4 @@
 {
-    rho1 = thermo1.rho();
-    rho2 = thermo2.rho();
-
     surfaceScalarField alpha1f(fvc::interpolate(alpha1));
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
 
@@ -91,10 +88,7 @@
     tmp<fvScalarMatrix> pEqnComp1;
     tmp<fvScalarMatrix> pEqnComp2;
 
-    //if (transonic)
-    //{
-    //}
-    //else
+    if (pimple.transonic())
     {
         surfaceScalarField phid1
         (
@@ -107,17 +101,42 @@
             fvc::interpolate(psi2)*phi2
         );
 
+        pEqnComp1 =
+            fvc::ddt(rho1)
+          + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1)
+          + correction
+            (
+                psi1*fvm::ddt(p)
+              + fvm::div(phid1, p) - fvm::Sp(fvc::div(phid1), p)
+            );
+        deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
+        pEqnComp1().relax();
+
+        pEqnComp2 =
+            fvc::ddt(rho2)
+          + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2)
+          + correction
+            (
+                psi2*fvm::ddt(p)
+              + fvm::div(phid2, p) - fvm::Sp(fvc::div(phid2), p)
+            );
+        deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
+        pEqnComp2().relax();
+    }
+    else
+    {
         pEqnComp1 =
             fvc::ddt(rho1) + psi1*correction(fvm::ddt(p))
-          + fvc::div(phid1, p)
-          - fvc::Sp(fvc::div(phid1), p);
+          + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1);
 
         pEqnComp2 =
             fvc::ddt(rho2) + psi2*correction(fvm::ddt(p))
-          + fvc::div(phid2, p)
-          - fvc::Sp(fvc::div(phid2), p);
+          + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2);
     }
 
+    // Cache p prior to solve for density update
+    volScalarField p_0(p);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqnIncomp
@@ -182,10 +201,9 @@
 
     p = max(p, pMin);
 
-    thermo1.correct();
-    thermo2.correct();
-    rho1 = thermo1.rho();
-    rho2 = thermo2.rho();
+    // Update densities from change in p
+    rho1 += psi1*(p - p_0);
+    rho2 += psi2*(p - p_0);
 
     K1 = 0.5*magSqr(U1);
     K2 = 0.5*magSqr(U2);
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H
index 0e373e1f404f7f747795202797049a0b212f367b..c4cdbc044bc772744e1769cce6ea5c4873de93c8 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H
@@ -53,6 +53,8 @@
             fvc::makeRelative(phi, U);
         }
     }
-
-    #include "continuityErrs.H"
 }
+
+phi.oldTime() = phi;
+
+#include "continuityErrs.H"
diff --git a/src/OSspecific/POSIX/signals/sigFpe.C b/src/OSspecific/POSIX/signals/sigFpe.C
index a275fd7a3ad34cb7951bdc76d110451420dea5c4..501ff95c4e4063fae3fdf38ba1b2491b2b3dbc87 100644
--- a/src/OSspecific/POSIX/signals/sigFpe.C
+++ b/src/OSspecific/POSIX/signals/sigFpe.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -263,8 +263,8 @@ void Foam::sigFpe::set(const bool verbose)
         {
             if (supported)
             {
-            Info<< "SetNaN : Initialising allocated memory to NaN"
-                << " (FOAM_SETNAN)." << endl;
+                Info<< "SetNaN : Initialising allocated memory to NaN"
+                    << " (FOAM_SETNAN)." << endl;
             }
             else
             {
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index dc8e946baaa521831ae24613b400dc5f9d249d78..033dd4f62c9fa5d76b5aef01198c171fc6b0e385 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -367,6 +367,8 @@ $(cellShape)/cellShapeIOList.C
 
 meshes/Identifiers/patch/patchIdentifier.C
 
+meshes/MeshObject/meshObject.C
+
 polyMesh = meshes/polyMesh
 
 polyPatches = $(polyMesh)/polyPatches
diff --git a/src/OpenFOAM/containers/Lists/PtrList/PtrListI.H b/src/OpenFOAM/containers/Lists/PtrList/PtrListI.H
index 3f29e9238f3cfa971f5711667501d76d1e15f45e..979e862d4a7c59e8eeb26e69d741bcbf750cd19e 100644
--- a/src/OpenFOAM/containers/Lists/PtrList/PtrListI.H
+++ b/src/OpenFOAM/containers/Lists/PtrList/PtrListI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -132,7 +132,9 @@ inline const T& Foam::PtrList<T>::operator[](const label i) const
     if (!ptrs_[i])
     {
         FatalErrorIn("PtrList::operator[] const")
-            << "hanging pointer, cannot dereference"
+            << "hanging pointer at index " << i
+            << " (size " << size()
+            << "), cannot dereference"
             << abort(FatalError);
     }
 
@@ -146,7 +148,9 @@ inline T& Foam::PtrList<T>::operator[](const label i)
     if (!ptrs_[i])
     {
         FatalErrorIn("PtrList::operator[]")
-            << "hanging pointer, cannot dereference"
+            << "hanging pointer at index " << i
+            << " (size " << size()
+            << "), cannot dereference"
             << abort(FatalError);
     }
 
diff --git a/src/OpenFOAM/containers/Lists/UPtrList/UPtrListI.H b/src/OpenFOAM/containers/Lists/UPtrList/UPtrListI.H
index dceeb1f279689244ad354d32499c951b5f43450e..38bb43093148519ad6353b9fd4fd423267fa3ff5 100644
--- a/src/OpenFOAM/containers/Lists/UPtrList/UPtrListI.H
+++ b/src/OpenFOAM/containers/Lists/UPtrList/UPtrListI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,7 +107,9 @@ inline const T& Foam::UPtrList<T>::operator[](const label i) const
     if (!ptrs_[i])
     {
         FatalErrorIn("UPtrList::operator[] const")
-            << "hanging pointer, cannot dereference"
+            << "hanging pointer at index " << i
+            << " (size " << size()
+            << "), cannot dereference"
             << abort(FatalError);
     }
 
@@ -121,7 +123,9 @@ inline T& Foam::UPtrList<T>::operator[](const label i)
     if (!ptrs_[i])
     {
         FatalErrorIn("UPtrList::operator[]")
-            << "hanging pointer, cannot dereference"
+            << "hanging pointer at index " << i
+            << " (size " << size()
+            << "), cannot dereference"
             << abort(FatalError);
     }
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C
index c3d68e28374dcb65032f72bf0befa9bcb8cdcb25..8a647b6db0392dbfa6e450781600ce2cbc209fce 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "processorPointPatchField.H"
-//#include "transformField.H"
 #include "processorPolyPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -92,73 +91,6 @@ processorPointPatchField<Type>::~processorPointPatchField()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-void processorPointPatchField<Type>::initSwapAddSeparated
-(
-    const Pstream::commsTypes commsType,
-    Field<Type>& pField
-)
-const
-{
-//    if (Pstream::parRun())
-//    {
-//        // Get internal field into correct order for opposite side
-//        Field<Type> pf
-//        (
-//            this->patchInternalField
-//            (
-//                pField,
-//                procPatch_.reverseMeshPoints()
-//            )
-//        );
-//
-//        OPstream::write
-//        (
-//            commsType,
-//            procPatch_.neighbProcNo(),
-//            reinterpret_cast<const char*>(pf.begin()),
-//            pf.byteSize(),
-//            procPatch_.tag()
-//        );
-//    }
-}
-
-
-template<class Type>
-void processorPointPatchField<Type>::swapAddSeparated
-(
-    const Pstream::commsTypes commsType,
-    Field<Type>& pField
-) const
-{
-//    if (Pstream::parRun())
-//    {
-//        Field<Type> pnf(this->size());
-//
-//        IPstream::read
-//        (
-//            commsType,
-//            procPatch_.neighbProcNo(),
-//            reinterpret_cast<char*>(pnf.begin()),
-//            pnf.byteSize(),
-//            procPatch_.tag()
-//        );
-//
-//        if (doTransform())
-//        {
-//            const processorPolyPatch& ppp = procPatch_.procPolyPatch();
-//            const tensor& forwardT = ppp.forwardT();
-//
-//            transform(pnf, forwardT, pnf);
-//        }
-//
-//        addToInternalField(pField, pnf, procPatch_.separatedPoints());
-//    }
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.H b/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.H
index 6582b62861afe20db4a12ccc3723de4c032370a2..5dfceda2d118ee6195f58f816a4affb2ed292144 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.H
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,7 +57,6 @@ class processorPointPatchField
         //- Local reference to processor patch
         const processorPointPatch& procPatch_;
 
-
 public:
 
     //- Runtime type information
@@ -176,19 +175,13 @@ public:
             )
             {}
 
-            //- Initialise swap of non-collocated patch point values
-            virtual void initSwapAddSeparated
-            (
-                const Pstream::commsTypes commsType,
-                Field<Type>&
-            ) const;
-
-            //- Complete swap of patch point values and add to local values
+            //- Assume processor patch always collocated
             virtual void swapAddSeparated
             (
                 const Pstream::commsTypes commsType,
                 Field<Type>&
-            ) const;
+            ) const
+            {}
 };
 
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/processorCyclic/processorCyclicPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/processorCyclic/processorCyclicPointPatchField.C
index 4e28c3466df64c53627cb2d7a6a14e51e387fe7a..6bea99a1bc998c4598358bb620f4b3c87b9e66c9 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/processorCyclic/processorCyclicPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/processorCyclic/processorCyclicPointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,8 @@ Foam::processorCyclicPointPatchField<Type>::processorCyclicPointPatchField
 )
 :
     coupledPointPatchField<Type>(p, iF),
-    procPatch_(refCast<const processorCyclicPointPatch>(p))
+    procPatch_(refCast<const processorCyclicPointPatch>(p)),
+    receiveBuf_(0)
 {}
 
 
@@ -51,7 +52,8 @@ Foam::processorCyclicPointPatchField<Type>::processorCyclicPointPatchField
 )
 :
     coupledPointPatchField<Type>(p, iF, dict),
-    procPatch_(refCast<const processorCyclicPointPatch>(p))
+    procPatch_(refCast<const processorCyclicPointPatch>(p)),
+    receiveBuf_(0)
 {}
 
 
@@ -65,7 +67,8 @@ Foam::processorCyclicPointPatchField<Type>::processorCyclicPointPatchField
 )
 :
     coupledPointPatchField<Type>(ptf, p, iF, mapper),
-    procPatch_(refCast<const processorCyclicPointPatch>(ptf.patch()))
+    procPatch_(refCast<const processorCyclicPointPatch>(ptf.patch())),
+    receiveBuf_(0)
 {}
 
 
@@ -77,7 +80,8 @@ Foam::processorCyclicPointPatchField<Type>::processorCyclicPointPatchField
 )
 :
     coupledPointPatchField<Type>(ptf, iF),
-    procPatch_(refCast<const processorCyclicPointPatch>(ptf.patch()))
+    procPatch_(refCast<const processorCyclicPointPatch>(ptf.patch())),
+    receiveBuf_(0)
 {}
 
 
@@ -109,6 +113,18 @@ void Foam::processorCyclicPointPatchField<Type>::initSwapAddSeparated
             )
         );
 
+        if (commsType == Pstream::nonBlocking)
+        {
+            receiveBuf_.setSize(pf.size());
+            IPstream::read
+            (
+                commsType,
+                procPatch_.neighbProcNo(),
+                reinterpret_cast<char*>(receiveBuf_.begin()),
+                receiveBuf_.byteSize(),
+                procPatch_.tag()
+            );
+        }
         OPstream::write
         (
             commsType,
@@ -130,16 +146,19 @@ void Foam::processorCyclicPointPatchField<Type>::swapAddSeparated
 {
     if (Pstream::parRun())
     {
-        Field<Type> pnf(this->size());
-
-        IPstream::read
-        (
-            commsType,
-            procPatch_.neighbProcNo(),
-            reinterpret_cast<char*>(pnf.begin()),
-            pnf.byteSize(),
-            procPatch_.tag()
-        );
+        // If nonblocking data has already been received into receiveBuf_
+        if (commsType != Pstream::nonBlocking)
+        {
+            receiveBuf_.setSize(this->size());
+            IPstream::read
+            (
+                commsType,
+                procPatch_.neighbProcNo(),
+                reinterpret_cast<char*>(receiveBuf_.begin()),
+                receiveBuf_.byteSize(),
+                procPatch_.tag()
+            );
+        }
 
         if (doTransform())
         {
@@ -147,11 +166,11 @@ void Foam::processorCyclicPointPatchField<Type>::swapAddSeparated
                 procPatch_.procCyclicPolyPatch();
             const tensor& forwardT = ppp.forwardT()[0];
 
-            transform(pnf, forwardT, pnf);
+            transform(receiveBuf_, forwardT, receiveBuf_);
         }
 
         // All points are separated
-        this->addToInternalField(pField, pnf);
+        this->addToInternalField(pField, receiveBuf_);
     }
 }
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/processorCyclic/processorCyclicPointPatchField.H b/src/OpenFOAM/fields/pointPatchFields/constraint/processorCyclic/processorCyclicPointPatchField.H
index c8575240603f5daeb6ba12b4fd5477e0cd53df66..6aa358b868b33f718f8fd717a3aab24141aa2489 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/processorCyclic/processorCyclicPointPatchField.H
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/processorCyclic/processorCyclicPointPatchField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,6 +57,9 @@ class processorCyclicPointPatchField
         //- Local reference to processor patch
         const processorCyclicPointPatch& procPatch_;
 
+        //- Receive buffer for non-blocking communication
+        mutable Field<Type> receiveBuf_;
+
 
 public:
 
diff --git a/src/OpenFOAM/meshes/MeshObject/MeshObject.C b/src/OpenFOAM/meshes/MeshObject/MeshObject.C
index 8f25049242ab9fcb1ce33f5458ede1f52c20ef0b..01e86da379045050f220dc211d2c4e5ba2a7beb5 100644
--- a/src/OpenFOAM/meshes/MeshObject/MeshObject.C
+++ b/src/OpenFOAM/meshes/MeshObject/MeshObject.C
@@ -25,6 +25,7 @@ License
 
 #include "MeshObject.H"
 #include "objectRegistry.H"
+#include "IOstreams.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -59,6 +60,11 @@ const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
     }
     else
     {
+        if (meshObject::debug)
+        {
+            Pout<< "MeshObject::New(const Mesh&) : constructing new "
+                << Type::typeName << endl;
+        }
         return regIOobject::store(new Type(mesh));
     }
 }
@@ -87,6 +93,11 @@ const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
     }
     else
     {
+        if (meshObject::debug)
+        {
+            Pout<< "MeshObject::New(const Mesh&) : constructing new "
+                << Type::typeName << endl;
+        }
         return regIOobject::store(new Type(mesh, d));
     }
 }
@@ -116,6 +127,11 @@ const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
     }
     else
     {
+        if (meshObject::debug)
+        {
+            Pout<< "MeshObject(const Mesh&) : constructing new "
+                << Type::typeName << endl;
+        }
         return regIOobject::store(new Type(mesh, d1, d2));
     }
 }
@@ -146,6 +162,11 @@ const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
     }
     else
     {
+        if (meshObject::debug)
+        {
+            Pout<< "MeshObject(const Mesh&) : constructing new "
+                << Type::typeName << endl;
+        }
         return regIOobject::store(new Type(mesh, d1, d2, d3));
     }
 }
@@ -177,6 +198,11 @@ const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
     }
     else
     {
+        if (meshObject::debug)
+        {
+            Pout<< "MeshObject(const Mesh&) : constructing new "
+                << Type::typeName << endl;
+        }
         return regIOobject::store(new Type(mesh, d1, d2, d3, d4));
     }
 }
@@ -195,6 +221,12 @@ bool Foam::MeshObject<Mesh, MeshObjectType, Type>::Delete(const Mesh& mesh)
         )
     )
     {
+        if (meshObject::debug)
+        {
+            Pout<< "MeshObject::Delete(const Mesh&) : deleting "
+                << Type::typeName << endl;
+        }
+
         return mesh.thisDb().checkOut
         (
             const_cast<Type&>
@@ -237,10 +269,22 @@ void Foam::meshObject::movePoints(objectRegistry& obr)
     {
         if (isA<MoveableMeshObject<Mesh> >(*iter()))
         {
+            if (meshObject::debug)
+            {
+                Pout<< "meshObject::movePoints(objectRegistry&) :"
+                    << " movePoints on "
+                    << iter()->name() << endl;
+            }
             dynamic_cast<MoveableMeshObject<Mesh>*>(iter())->movePoints();
         }
         else
         {
+            if (meshObject::debug)
+            {
+                Pout<< "meshObject::movePoints(objectRegistry&) :"
+                    << " destroying "
+                    << iter()->name() << endl;
+            }
             obr.checkOut(*iter());
         }
     }
@@ -264,10 +308,21 @@ void Foam::meshObject::updateMesh(objectRegistry& obr, const mapPolyMesh& mpm)
     {
         if (isA<UpdateableMeshObject<Mesh> >(*iter()))
         {
+            if (meshObject::debug)
+            {
+                Pout<< "meshObject::updateMesh(objectRegistry&) :"
+                    << " updateMesh on "
+                    << iter()->name() << endl;
+            }
             dynamic_cast<UpdateableMeshObject<Mesh>*>(iter())->updateMesh(mpm);
         }
         else
         {
+            if (meshObject::debug)
+            {
+                Pout<< "meshObject::updateMesh(objectRegistry&) : destroying "
+                    << iter()->name() << endl;
+            }
             obr.checkOut(*iter());
         }
     }
@@ -284,6 +339,11 @@ void Foam::meshObject::clear(objectRegistry& obr)
 
     forAllIter(typename HashTable<MeshObjectType<Mesh>*>, meshObjects, iter)
     {
+        if (meshObject::debug)
+        {
+            Pout<< "meshObject::clear(objectRegistry&) : destroying "
+                << iter()->name() << endl;
+        }
         obr.checkOut(*iter());
     }
 }
diff --git a/src/OpenFOAM/meshes/MeshObject/MeshObject.H b/src/OpenFOAM/meshes/MeshObject/MeshObject.H
index 8efa89f944e9ac5e3ddd5aa6a8ac292db0aa3e18..b95c631e52fe1147c297d3d2645ad93130fe6d2a 100644
--- a/src/OpenFOAM/meshes/MeshObject/MeshObject.H
+++ b/src/OpenFOAM/meshes/MeshObject/MeshObject.H
@@ -164,20 +164,12 @@ class meshObject
 {
 public:
 
+    // Declare name of the class and its debug switch
+    ClassName("meshObject");
+
     // Constructors
 
-        meshObject(const word& typeName, const objectRegistry& obr)
-        :
-            regIOobject
-            (
-                IOobject
-                (
-                    typeName,
-                    obr.instance(),
-                    obr
-                )
-            )
-        {}
+        meshObject(const word& typeName, const objectRegistry& obr);
 
 
     // Static member functions
diff --git a/src/OpenFOAM/meshes/MeshObject/meshObject.C b/src/OpenFOAM/meshes/MeshObject/meshObject.C
new file mode 100644
index 0000000000000000000000000000000000000000..7130f38b5dfd738ff2d9ef9b4c6f2316dc24badd
--- /dev/null
+++ b/src/OpenFOAM/meshes/MeshObject/meshObject.C
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MeshObject.H"
+
+/* * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * */
+
+namespace Foam
+{
+    defineTypeNameAndDebug(meshObject, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::meshObject::meshObject(const word& typeName, const objectRegistry& obr)
+:
+    regIOobject
+    (
+        IOobject
+        (
+            typeName,
+            obr.instance(),
+            obr
+        )
+    )
+{}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C
index 8988a06eb211dad9cf592dddbdec24def09fe970..7e45266d241298f134d0232c163651160fafadda 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.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
@@ -37,10 +37,7 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
 (
@@ -66,20 +63,10 @@ Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
         neighbProcNo,
         transform
     ),
-    tag_
-    (
-        Pstream::nProcs()*max(myProcNo, neighbProcNo)
-      + min(myProcNo, neighbProcNo)
-    ),
     referPatchName_(referPatchName),
+    tag_(-1),
     referPatchID_(-1)
-{
-    if (debug)
-    {
-        Pout<< "processorCyclicPolyPatch " << name << " uses tag " << tag_
-            << endl;
-    }
-}
+{}
 
 
 Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
@@ -92,20 +79,10 @@ Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
 )
 :
     processorPolyPatch(name, dict, index, bm, patchType),
-    tag_
-    (
-        Pstream::nProcs()*max(myProcNo(), neighbProcNo())
-      + min(myProcNo(), neighbProcNo())
-    ),
     referPatchName_(dict.lookup("referPatch")),
+    tag_(dict.lookupOrDefault<int>("tag", -1)),
     referPatchID_(-1)
-{
-    if (debug)
-    {
-        Pout<< "processorCyclicPolyPatch " << name << " uses tag " << tag_
-            << endl;
-    }
-}
+{}
 
 
 Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
@@ -115,8 +92,8 @@ Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
 )
 :
     processorPolyPatch(pp, bm),
-    tag_(pp.tag_),
     referPatchName_(pp.referPatchName()),
+    tag_(pp.tag()),
     referPatchID_(-1)
 {}
 
@@ -132,8 +109,8 @@ Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
 )
 :
     processorPolyPatch(pp, bm, index, newSize, newStart),
-    tag_(pp.tag_),
     referPatchName_(referPatchName),
+    tag_(-1),
     referPatchID_(-1)
 {}
 
@@ -148,8 +125,8 @@ Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
 )
 :
     processorPolyPatch(pp, bm, index, mapAddressing, newStart),
-    tag_(pp.tag_),
     referPatchName_(pp.referPatchName()),
+    tag_(-1),
     referPatchID_(-1)
 {}
 
@@ -162,6 +139,45 @@ Foam::processorCyclicPolyPatch::~processorCyclicPolyPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+int Foam::processorCyclicPolyPatch::tag() const
+{
+    if (tag_ == -1)
+    {
+        // Get unique tag to use for all comms. Make sure that both sides
+        // use the same tag
+        const cyclicPolyPatch& cycPatch = refCast<const cyclicPolyPatch>
+        (
+            referPatch()
+        );
+
+        if (cycPatch.owner())
+        {
+            tag_ = Hash<word>()(cycPatch.name());
+        }
+        else
+        {
+            tag_ = Hash<word>()(cycPatch.neighbPatch().name());
+        }
+
+        if (tag_ == Pstream::msgType() || tag_ == -1)
+        {
+            FatalErrorIn("processorCyclicPolyPatch::tag() const")
+                << "Tag calculated from cyclic patch name " << tag_
+                << " is the same as the current message type "
+                << Pstream::msgType() << " or -1" << nl
+                << "Please set a non-conflicting, unique, tag by hand"
+                << " using the 'tag' entry"
+                << exit(FatalError);
+        }
+        if (debug)
+        {
+            Pout<< "processorCyclicPolyPatch " << name() << " uses tag " << tag_
+                << endl;
+        }
+    }
+    return tag_;
+}
+
 
 void Foam::processorCyclicPolyPatch::initGeometry(PstreamBuffers& pBufs)
 {
@@ -283,6 +299,11 @@ void Foam::processorCyclicPolyPatch::write(Ostream& os) const
     processorPolyPatch::write(os);
     os.writeKeyword("referPatch") << referPatchName_
         << token::END_STATEMENT << nl;
+    if (tag_ != -1)
+    {
+        os.writeKeyword("tag") << tag_
+            << token::END_STATEMENT << nl;
+    }
 }
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.H
index fa59414586b676ac17899a43c6a822d384b009aa..0b9d4f70d94439fcf120b74ea404f75105a5aa53 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.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
@@ -54,12 +54,12 @@ class processorCyclicPolyPatch
 {
     // Private data
 
-        //- Message tag to use for communication
-        const int tag_;
-
         //- Name of originating patch
         const word referPatchName_;
 
+        //- Message tag to use for communication
+        mutable int tag_;
+
         //- Index of originating patch
         mutable label referPatchID_;
 
@@ -264,10 +264,7 @@ public:
         }
 
         //- Return message tag to use for communication
-        virtual int tag() const
-        {
-            return tag_;
-        }
+        virtual int tag() const;
 
         //- Does this side own the patch ?
         virtual bool owner() const
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C
index bd15d4764ad68e6b1de1b23dcb2558de80e89348..5d359bc691550b2a9c65b443acf2d5bbf4937571 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,12 +42,9 @@ mappedFieldFvPatchField<Type>::mappedFieldFvPatchField
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    mappedPatchBase(p.patch()),
     fixedValueFvPatchField<Type>(p, iF),
-    fieldName_(iF.name()),
-    setAverage_(false),
-    average_(pTraits<Type>::zero),
-    interpolationScheme_(interpolationCell<Type>::typeName)
+    mappedPatchBase(p.patch()),
+    mappedPatchFieldBase<Type>(*this, *this)
 {}
 
 
@@ -60,12 +57,9 @@ mappedFieldFvPatchField<Type>::mappedFieldFvPatchField
     const fvPatchFieldMapper& mapper
 )
 :
-    mappedPatchBase(p.patch(), ptf),
     fixedValueFvPatchField<Type>(ptf, p, iF, mapper),
-    fieldName_(ptf.fieldName_),
-    setAverage_(ptf.setAverage_),
-    average_(ptf.average_),
-    interpolationScheme_(ptf.interpolationScheme_)
+    mappedPatchBase(p.patch(), ptf),
+    mappedPatchFieldBase<Type>(*this, *this, ptf)
 {}
 
 
@@ -77,18 +71,10 @@ mappedFieldFvPatchField<Type>::mappedFieldFvPatchField
     const dictionary& dict
 )
 :
-    mappedPatchBase(p.patch(), dict),
     fixedValueFvPatchField<Type>(p, iF, dict),
-    fieldName_(dict.template lookupOrDefault<word>("fieldName", iF.name())),
-    setAverage_(readBool(dict.lookup("setAverage"))),
-    average_(pTraits<Type>(dict.lookup("average"))),
-    interpolationScheme_(interpolationCell<Type>::typeName)
-{
-    if (mode() == mappedPatchBase::NEARESTCELL)
-    {
-        dict.lookup("interpolationScheme") >> interpolationScheme_;
-    }
-}
+    mappedPatchBase(p.patch(), dict),
+    mappedPatchFieldBase<Type>(*this, *this, dict)
+{}
 
 
 template<class Type>
@@ -110,6 +96,7 @@ mappedFieldFvPatchField<Type>::mappedFieldFvPatchField
     const word& interpolationScheme
 )
 :
+    fixedValueFvPatchField<Type>(p, iF),
     mappedPatchBase
     (
         p.patch(),
@@ -118,11 +105,15 @@ mappedFieldFvPatchField<Type>::mappedFieldFvPatchField
         samplePatch,
         distance
     ),
-    fixedValueFvPatchField<Type>(p, iF),
-    fieldName_(fieldName),
-    setAverage_(setAverage),
-    average_(average),
-    interpolationScheme_(interpolationScheme)
+    mappedPatchFieldBase<Type>
+    (
+        *this,
+        *this,
+        fieldName,
+        setAverage,
+        average,
+        interpolationScheme
+    )
 {}
 
 
@@ -132,12 +123,9 @@ mappedFieldFvPatchField<Type>::mappedFieldFvPatchField
     const mappedFieldFvPatchField<Type>& ptf
 )
 :
-    mappedPatchBase(ptf.patch().patch(), ptf),
     fixedValueFvPatchField<Type>(ptf),
-    fieldName_(ptf.fieldName_),
-    setAverage_(ptf.setAverage_),
-    average_(ptf.average_),
-    interpolationScheme_(ptf.interpolationScheme_)
+    mappedPatchBase(ptf.patch().patch(), ptf),
+    mappedPatchFieldBase<Type>(ptf)
 {}
 
 
@@ -148,65 +136,14 @@ mappedFieldFvPatchField<Type>::mappedFieldFvPatchField
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    mappedPatchBase(ptf.patch().patch(), ptf),
     fixedValueFvPatchField<Type>(ptf, iF),
-    fieldName_(ptf.fieldName_),
-    setAverage_(ptf.setAverage_),
-    average_(ptf.average_),
-    interpolationScheme_(ptf.interpolationScheme_)
+    mappedPatchBase(ptf.patch().patch(), ptf),
+    mappedPatchFieldBase<Type>(*this, *this, ptf)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type>
-const GeometricField<Type, fvPatchField, volMesh>&
-mappedFieldFvPatchField<Type>::sampleField() const
-{
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    const fvMesh& nbrMesh = refCast<const fvMesh>(sampleMesh());
-
-    if (sameRegion())
-    {
-        if (fieldName_ == this->dimensionedInternalField().name())
-        {
-            // Optimisation: bypass field lookup
-            return
-                dynamic_cast<const fieldType&>
-                (
-                    this->dimensionedInternalField()
-                );
-        }
-        else
-        {
-            const fvMesh& thisMesh = this->patch().boundaryMesh().mesh();
-            return thisMesh.template lookupObject<fieldType>(fieldName_);
-        }
-    }
-    else
-    {
-        return nbrMesh.template lookupObject<fieldType>(fieldName_);
-    }
-}
-
-
-template<class Type>
-const interpolation<Type>&
-mappedFieldFvPatchField<Type>::interpolator() const
-{
-    if (!interpolator_.valid())
-    {
-        interpolator_ = interpolation<Type>::New
-        (
-            interpolationScheme_,
-            sampleField()
-        );
-    }
-    return interpolator_();
-}
-
-
 template<class Type>
 void mappedFieldFvPatchField<Type>::updateCoeffs()
 {
@@ -215,132 +152,7 @@ void mappedFieldFvPatchField<Type>::updateCoeffs()
         return;
     }
 
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    // Since we're inside initEvaluate/evaluate there might be processor
-    // comms underway. Change the tag we use.
-    int oldTag = UPstream::msgType();
-    UPstream::msgType() = oldTag + 1;
-
-    const fvMesh& thisMesh = this->patch().boundaryMesh().mesh();
-    const fvMesh& nbrMesh = refCast<const fvMesh>(sampleMesh());
-
-    // Result of obtaining remote values
-    Field<Type> newValues;
-
-    switch (mode())
-    {
-        case NEARESTCELL:
-        {
-            const mapDistribute& mapDist = this->mappedPatchBase::map();
-
-            if (interpolationScheme_ != interpolationCell<Type>::typeName)
-            {
-                // Need to do interpolation so need cells to sample
-
-                // Send back sample points to the processor that holds the cell
-                vectorField samples(samplePoints());
-                mapDist.reverseDistribute
-                (
-                    (sameRegion() ? thisMesh.nCells() : nbrMesh.nCells()),
-                    point::max,
-                    samples
-                );
-
-                const interpolation<Type>& interp = interpolator();
-
-                newValues.setSize(samples.size(), pTraits<Type>::max);
-                forAll(samples, cellI)
-                {
-                    if (samples[cellI] != point::max)
-                    {
-                        newValues[cellI] = interp.interpolate
-                        (
-                            samples[cellI],
-                            cellI
-                        );
-                    }
-                }
-            }
-            else
-            {
-                newValues = sampleField();
-            }
-
-            mapDist.distribute(newValues);
-
-            break;
-        }
-        case NEARESTPATCHFACE: case NEARESTPATCHFACEAMI:
-        {
-            const label nbrPatchID =
-                nbrMesh.boundaryMesh().findPatchID(samplePatch());
-            if (nbrPatchID < 0)
-            {
-                FatalErrorIn
-                (
-                    "void mappedFieldFvPatchField<Type>::updateCoeffs()"
-                )<< "Unable to find sample patch " << samplePatch()
-                 << " in region " << sampleRegion()
-                 << " for patch " << this->patch().name() << nl
-                 << abort(FatalError);
-            }
-
-            const fieldType& nbrField = sampleField();
-
-            newValues = nbrField.boundaryField()[nbrPatchID];
-            this->distribute(newValues);
-
-            break;
-        }
-        case NEARESTFACE:
-        {
-            Field<Type> allValues(nbrMesh.nFaces(), pTraits<Type>::zero);
-
-            const fieldType& nbrField = sampleField();
-
-            forAll(nbrField.boundaryField(), patchI)
-            {
-                const fvPatchField<Type>& pf =
-                    nbrField.boundaryField()[patchI];
-                label faceStart = pf.patch().patch().start();
-
-                forAll(pf, faceI)
-                {
-                    allValues[faceStart++] = pf[faceI];
-                }
-            }
-
-            this->distribute(allValues);
-            newValues.transfer(allValues);
-
-            break;
-        }
-        default:
-        {
-            FatalErrorIn("mappedFieldFvPatchField<Type>::updateCoeffs()")
-                << "Unknown sampling mode: " << mode()
-                << nl << abort(FatalError);
-        }
-    }
-
-    if (setAverage_)
-    {
-        Type averagePsi =
-            gSum(this->patch().magSf()*newValues)
-           /gSum(this->patch().magSf());
-
-        if (mag(averagePsi)/mag(average_) > 0.5)
-        {
-            newValues *= mag(average_)/mag(averagePsi);
-        }
-        else
-        {
-            newValues += (average_ - averagePsi);
-        }
-    }
-
-    this->operator==(newValues);
+    this->operator==(this->mappedField());
 
     if (debug)
     {
@@ -352,9 +164,6 @@ void mappedFieldFvPatchField<Type>::updateCoeffs()
             << endl;
     }
 
-    // Restore tag
-    UPstream::msgType() = oldTag;
-
     fixedValueFvPatchField<Type>::updateCoeffs();
 }
 
@@ -364,11 +173,7 @@ void mappedFieldFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
     mappedPatchBase::write(os);
-    os.writeKeyword("fieldName") << fieldName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl;
-    os.writeKeyword("average") << average_ << token::END_STATEMENT << nl;
-    os.writeKeyword("interpolationScheme") << interpolationScheme_
-        << token::END_STATEMENT << nl;
+    mappedPatchFieldBase<Type>::write(os);
     this->writeEntry("value", os);
 }
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.H
index 6dd9f3efb5c8579b8c9c966394d4aa34a5e8685a..983df12eb4aec8ecf35dceb4476bd3922c279613 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.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
@@ -75,6 +75,7 @@ SourceFiles
 #define mappedFieldFvPatchField_H
 
 #include "mappedPatchBase.H"
+#include "mappedPatchFieldBase.H"
 #include "fixedValueFvPatchFields.H"
 #include "interpolation.H"
 
@@ -90,37 +91,10 @@ namespace Foam
 template<class Type>
 class mappedFieldFvPatchField
 :
+    public fixedValueFvPatchField<Type>,
     public mappedPatchBase,
-    public fixedValueFvPatchField<Type>
+    public mappedPatchFieldBase<Type>
 {
-    // Private data
-
-        //- Name of field to sample - defaults to field associated with this
-        //  patchField if not specified
-        word fieldName_;
-
-        //- If true adjust the mapped field to maintain average value average_
-        const bool setAverage_;
-
-        //- Average value the mapped field is adjusted to maintain if
-        //  setAverage_ is set true
-        const Type average_;
-
-        //- Interpolation scheme to use for nearestCell mode
-        word interpolationScheme_;
-
-        //- Pointer to the cell interpolator
-        mutable autoPtr<interpolation<Type> > interpolator_;
-
-
-    // Private Member Functions
-
-        //- Field to sample. Either on my or nbr mesh
-        const GeometricField<Type, fvPatchField, volMesh>& sampleField() const;
-
-        //- Access the interpolation method
-        const interpolation<Type>& interpolator() const;
-
 
 public:
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.C
new file mode 100644
index 0000000000000000000000000000000000000000..291fe316dce70579cc6353b3380f1713888f2990
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.C
@@ -0,0 +1,340 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "mappedPatchFieldBase.H"
+#include "mappedPatchBase.H"
+#include "interpolationCell.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+mappedPatchFieldBase<Type>::mappedPatchFieldBase
+(
+    const mappedPatchBase& mapper,
+    const fvPatchField<Type>& patchField,
+    const word& fieldName,
+    const bool setAverage,
+    const Type average,
+    const word& interpolationScheme
+)
+:
+    mapper_(mapper),
+    patchField_(patchField),
+    fieldName_(fieldName),
+    setAverage_(setAverage),
+    average_(average),
+    interpolationScheme_(interpolationScheme)
+{}
+
+
+template<class Type>
+mappedPatchFieldBase<Type>::mappedPatchFieldBase
+(
+    const mappedPatchBase& mapper,
+    const fvPatchField<Type>& patchField,
+    const dictionary& dict
+)
+:
+    mapper_(mapper),
+    patchField_(patchField),
+    fieldName_
+    (
+        dict.template lookupOrDefault<word>
+        (
+            "fieldName",
+            patchField_.dimensionedInternalField().name()
+        )
+    ),
+    setAverage_(readBool(dict.lookup("setAverage"))),
+    average_(pTraits<Type>(dict.lookup("average"))),
+    interpolationScheme_(interpolationCell<Type>::typeName)
+{
+    if (mapper_.mode() == mappedPatchBase::NEARESTCELL)
+    {
+        dict.lookup("interpolationScheme") >> interpolationScheme_;
+    }
+}
+
+
+template<class Type>
+mappedPatchFieldBase<Type>::mappedPatchFieldBase
+(
+    const mappedPatchBase& mapper,
+    const fvPatchField<Type>& patchField
+)
+:
+    mapper_(mapper),
+    patchField_(patchField),
+    fieldName_(patchField_.dimensionedInternalField().name()),
+    setAverage_(false),
+    average_(pTraits<Type>::zero),
+    interpolationScheme_(interpolationCell<Type>::typeName)
+{}
+
+
+template<class Type>
+mappedPatchFieldBase<Type>::mappedPatchFieldBase
+(
+    const mappedPatchFieldBase<Type>& mapper
+)
+:
+    mapper_(mapper.mapper_),
+    patchField_(mapper.patchField_),
+    fieldName_(mapper.fieldName_),
+    setAverage_(mapper.setAverage_),
+    average_(mapper.average_),
+    interpolationScheme_(mapper.interpolationScheme_)
+{}
+
+
+template<class Type>
+mappedPatchFieldBase<Type>::mappedPatchFieldBase
+(
+    const mappedPatchBase& mapper,
+    const fvPatchField<Type>& patchField,
+    const mappedPatchFieldBase<Type>& base
+)
+:
+    mapper_(mapper),
+    patchField_(patchField),
+    fieldName_(base.fieldName_),
+    setAverage_(base.setAverage_),
+    average_(base.average_),
+    interpolationScheme_(base.interpolationScheme_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+const GeometricField<Type, fvPatchField, volMesh>&
+mappedPatchFieldBase<Type>::sampleField() const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    const fvMesh& nbrMesh = refCast<const fvMesh>(mapper_.sampleMesh());
+
+    if (mapper_.sameRegion())
+    {
+        if (fieldName_ == patchField_.dimensionedInternalField().name())
+        {
+            // Optimisation: bypass field lookup
+            return
+                dynamic_cast<const fieldType&>
+                (
+                    patchField_.dimensionedInternalField()
+                );
+        }
+        else
+        {
+            const fvMesh& thisMesh = patchField_.patch().boundaryMesh().mesh();
+            return thisMesh.template lookupObject<fieldType>(fieldName_);
+        }
+    }
+    else
+    {
+        return nbrMesh.template lookupObject<fieldType>(fieldName_);
+    }
+}
+
+
+template<class Type>
+const interpolation<Type>& mappedPatchFieldBase<Type>::interpolator() const
+{
+    if (!interpolator_.valid())
+    {
+        interpolator_ = interpolation<Type>::New
+        (
+            interpolationScheme_,
+            sampleField()
+        );
+    }
+    return interpolator_();
+}
+
+
+template<class Type>
+tmp<Field<Type> > mappedPatchFieldBase<Type>::mappedField() const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    // Since we're inside initEvaluate/evaluate there might be processor
+    // comms underway. Change the tag we use.
+    int oldTag = UPstream::msgType();
+    UPstream::msgType() = oldTag + 1;
+
+    const fvMesh& thisMesh = patchField_.patch().boundaryMesh().mesh();
+    const fvMesh& nbrMesh = refCast<const fvMesh>(mapper_.sampleMesh());
+
+    // Result of obtaining remote values
+    tmp<Field<Type> > tnewValues(new Field<Type>(0));
+    Field<Type>& newValues = tnewValues();
+
+    switch (mapper_.mode())
+    {
+        case mappedPatchBase::NEARESTCELL:
+        {
+            const mapDistribute& distMap = mapper_.map();
+
+            if (interpolationScheme_ != interpolationCell<Type>::typeName)
+            {
+                // Send back sample points to the processor that holds the cell
+                vectorField samples(mapper_.samplePoints());
+                distMap.reverseDistribute
+                (
+                    (
+                        mapper_.sameRegion()
+                      ? thisMesh.nCells()
+                      : nbrMesh.nCells()
+                    ),
+                    point::max,
+                    samples
+                );
+
+                const interpolation<Type>& interp = interpolator();
+
+                newValues.setSize(samples.size(), pTraits<Type>::max);
+                forAll(samples, cellI)
+                {
+                    if (samples[cellI] != point::max)
+                    {
+                        newValues[cellI] = interp.interpolate
+                        (
+                            samples[cellI],
+                            cellI
+                        );
+                    }
+                }
+            }
+            else
+            {
+                newValues = sampleField();
+            }
+
+            distMap.distribute(newValues);
+
+            break;
+        }
+        case mappedPatchBase::NEARESTPATCHFACE:
+        case mappedPatchBase::NEARESTPATCHFACEAMI:
+        {
+            const label nbrPatchID =
+                nbrMesh.boundaryMesh().findPatchID(mapper_.samplePatch());
+
+            if (nbrPatchID < 0)
+            {
+                FatalErrorIn
+                (
+                    "void mappedPatchFieldBase<Type>::updateCoeffs()"
+                )<< "Unable to find sample patch " << mapper_.samplePatch()
+                 << " in region " << mapper_.sampleRegion()
+                 << " for patch " << patchField_.patch().name() << nl
+                 << abort(FatalError);
+            }
+
+            const fieldType& nbrField = sampleField();
+
+            newValues = nbrField.boundaryField()[nbrPatchID];
+            mapper_.distribute(newValues);
+
+            break;
+        }
+        case mappedPatchBase::NEARESTFACE:
+        {
+            Field<Type> allValues(nbrMesh.nFaces(), pTraits<Type>::zero);
+
+            const fieldType& nbrField = sampleField();
+
+            forAll(nbrField.boundaryField(), patchI)
+            {
+                const fvPatchField<Type>& pf =
+                    nbrField.boundaryField()[patchI];
+                label faceStart = pf.patch().start();
+
+                forAll(pf, faceI)
+                {
+                    allValues[faceStart++] = pf[faceI];
+                }
+            }
+
+            mapper_.distribute(allValues);
+            newValues.transfer(allValues);
+
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "mappedPatchFieldBase<Type>::updateCoeffs()"
+            )<< "Unknown sampling mode: " << mapper_.mode()
+             << nl << abort(FatalError);
+        }
+    }
+
+    if (setAverage_)
+    {
+        Type averagePsi =
+            gSum(patchField_.patch().magSf()*newValues)
+           /gSum(patchField_.patch().magSf());
+
+        if (mag(averagePsi)/mag(average_) > 0.5)
+        {
+            newValues *= mag(average_)/mag(averagePsi);
+        }
+        else
+        {
+            newValues += (average_ - averagePsi);
+        }
+    }
+
+    // Restore tag
+    UPstream::msgType() = oldTag;
+
+    return tnewValues;
+}
+
+
+template<class Type>
+void mappedPatchFieldBase<Type>::write(Ostream& os) const
+{
+    os.writeKeyword("fieldName") << fieldName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl;
+    os.writeKeyword("average") << average_ << token::END_STATEMENT << nl;
+    os.writeKeyword("interpolationScheme") << interpolationScheme_
+        << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.H b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.H
new file mode 100644
index 0000000000000000000000000000000000000000..50728410949b23ae732b5110718dfc7a3b2642c3
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.H
@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::mappedPatchFieldBase
+
+Description
+    Functionality for sampling fields using mappedPatchBase.
+
+SourceFiles
+    mappedPatchFieldBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mappedPatchFieldBase_H
+#define mappedPatchFieldBase_H
+
+#include "fixedValueFvPatchFields.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class mappedPatchBase;
+template<class> class interpolation;
+
+/*---------------------------------------------------------------------------*\
+                  Class mappedPatchFieldBase Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class mappedPatchFieldBase
+{
+
+protected:
+
+    // Protected data
+
+        //- Mapping engine
+        const mappedPatchBase& mapper_;
+
+        //- Underlying patch field
+        const fvPatchField<Type>& patchField_;
+
+        //- Name of field to sample
+        word fieldName_;
+
+        //- If true adjust the mapped field to maintain average value average_
+        const bool setAverage_;
+
+        //- Average value the mapped field is adjusted to maintain if
+        //  setAverage_ is set true
+        const Type average_;
+
+        //- Interpolation scheme to use for nearestcell mode
+        word interpolationScheme_;
+
+        mutable autoPtr<interpolation<Type> > interpolator_;
+
+
+    // Protected Member Functions
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        mappedPatchFieldBase
+        (
+            const mappedPatchBase& mapper,
+            const fvPatchField<Type>& patchField,
+            const word& fieldName,
+            const bool setAverage,
+            const Type average,
+            const word& interpolationScheme
+        );
+
+        //- Construct from dictionary
+        mappedPatchFieldBase
+        (
+            const mappedPatchBase& mapper,
+            const fvPatchField<Type>& patchField,
+            const dictionary& dict
+        );
+
+        //- Construct empty
+        mappedPatchFieldBase
+        (
+            const mappedPatchBase& mapper,
+            const fvPatchField<Type>& patchField
+        );
+
+        //- Construct copy
+        mappedPatchFieldBase
+        (
+            const mappedPatchFieldBase<Type>& mapper
+        );
+
+        //- Construct copy, resetting patch and field
+        mappedPatchFieldBase
+        (
+            const mappedPatchBase& mapper,
+            const fvPatchField<Type>& patchField,
+            const mappedPatchFieldBase<Type>& base
+        );
+
+
+    //- Destructor
+    virtual ~mappedPatchFieldBase<Type>()
+    {}
+
+
+    // Member functions
+
+        //- Field to sample. Either on my or nbr mesh
+        const GeometricField<Type, fvPatchField, volMesh>& sampleField() const;
+
+        //- Access the interpolation method
+        const interpolation<Type>& interpolator() const;
+
+        //- Map sampleField onto *this patch
+        virtual tmp<Field<Type> > mappedField() const;
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "mappedPatchFieldBase.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.C
index 77cc8d61a9d8dad64d7cf0935df71d716bf1d167..851344194b554ac57457754983d74cddc289eb59 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.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,7 +26,6 @@ License
 #include "mappedFixedValueFvPatchField.H"
 #include "mappedPatchBase.H"
 #include "volFields.H"
-#include "interpolationCell.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,10 +42,7 @@ mappedFixedValueFvPatchField<Type>::mappedFixedValueFvPatchField
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    fieldName_(iF.name()),
-    setAverage_(false),
-    average_(pTraits<Type>::zero),
-    interpolationScheme_(interpolationCell<Type>::typeName)
+    mappedPatchFieldBase<Type>(this->mapper(p, iF), *this)
 {}
 
 
@@ -60,31 +56,8 @@ mappedFixedValueFvPatchField<Type>::mappedFixedValueFvPatchField
 )
 :
     fixedValueFvPatchField<Type>(ptf, p, iF, mapper),
-    fieldName_(ptf.fieldName_),
-    setAverage_(ptf.setAverage_),
-    average_(ptf.average_),
-    interpolationScheme_(ptf.interpolationScheme_)
-{
-    if (!isA<mappedPatchBase>(this->patch().patch()))
-    {
-        FatalErrorIn
-        (
-            "mappedFixedValueFvPatchField<Type>::"
-            "mappedFixedValueFvPatchField\n"
-            "(\n"
-            "    const mappedFixedValueFvPatchField<Type>&,\n"
-            "    const fvPatch&,\n"
-            "    const Field<Type>&,\n"
-            "    const fvPatchFieldMapper&\n"
-            ")\n"
-        )   << "\n    patch type '" << p.type()
-            << "' not type '" << mappedPatchBase::typeName << "'"
-            << "\n    for patch " << p.name()
-            << " of field " << this->dimensionedInternalField().name()
-            << " in file " << this->dimensionedInternalField().objectPath()
-            << exit(FatalError);
-    }
-}
+    mappedPatchFieldBase<Type>(this->mapper(p, iF), *this, ptf)
+{}
 
 
 template<class Type>
@@ -96,39 +69,8 @@ mappedFixedValueFvPatchField<Type>::mappedFixedValueFvPatchField
 )
 :
     fixedValueFvPatchField<Type>(p, iF, dict),
-    fieldName_(dict.lookupOrDefault<word>("fieldName", iF.name())),
-    setAverage_(readBool(dict.lookup("setAverage"))),
-    average_(pTraits<Type>(dict.lookup("average"))),
-    interpolationScheme_(interpolationCell<Type>::typeName)
-{
-    if (!isA<mappedPatchBase>(this->patch().patch()))
-    {
-        FatalErrorIn
-        (
-            "mappedFixedValueFvPatchField<Type>::"
-            "mappedFixedValueFvPatchField\n"
-            "(\n"
-            "    const fvPatch& p,\n"
-            "    const DimensionedField<Type, volMesh>& iF,\n"
-            "    const dictionary& dict\n"
-            ")\n"
-        )   << "\n    patch type '" << p.type()
-            << "' not type '" << mappedPatchBase::typeName << "'"
-            << "\n    for patch " << p.name()
-            << " of field " << this->dimensionedInternalField().name()
-            << " in file " << this->dimensionedInternalField().objectPath()
-            << exit(FatalError);
-    }
-
-    const mappedPatchBase& mpp = refCast<const mappedPatchBase>
-    (
-        mappedFixedValueFvPatchField<Type>::patch().patch()
-    );
-    if (mpp.mode() == mappedPatchBase::NEARESTCELL)
-    {
-        dict.lookup("interpolationScheme") >> interpolationScheme_;
-    }
-}
+    mappedPatchFieldBase<Type>(this->mapper(p, iF), *this, dict)
+{}
 
 
 template<class Type>
@@ -138,10 +80,7 @@ mappedFixedValueFvPatchField<Type>::mappedFixedValueFvPatchField
 )
 :
     fixedValueFvPatchField<Type>(ptf),
-    fieldName_(ptf.fieldName_),
-    setAverage_(ptf.setAverage_),
-    average_(ptf.average_),
-    interpolationScheme_(ptf.interpolationScheme_)
+    mappedPatchFieldBase<Type>(ptf)
 {}
 
 
@@ -153,64 +92,32 @@ mappedFixedValueFvPatchField<Type>::mappedFixedValueFvPatchField
 )
 :
     fixedValueFvPatchField<Type>(ptf, iF),
-    fieldName_(ptf.fieldName_),
-    setAverage_(ptf.setAverage_),
-    average_(ptf.average_),
-    interpolationScheme_(ptf.interpolationScheme_)
+    mappedPatchFieldBase<Type>(this->mapper(this->patch(), iF), *this, ptf)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-const GeometricField<Type, fvPatchField, volMesh>&
-mappedFixedValueFvPatchField<Type>::sampleField() const
-{
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    const mappedPatchBase& mpp = refCast<const mappedPatchBase>
-    (
-        mappedFixedValueFvPatchField<Type>::patch().patch()
-    );
-    const fvMesh& nbrMesh = refCast<const fvMesh>(mpp.sampleMesh());
-
-    if (mpp.sameRegion())
-    {
-        if (fieldName_ == this->dimensionedInternalField().name())
-        {
-            // Optimisation: bypass field lookup
-            return
-                dynamic_cast<const fieldType&>
-                (
-                    this->dimensionedInternalField()
-                );
-        }
-        else
-        {
-            const fvMesh& thisMesh = this->patch().boundaryMesh().mesh();
-            return thisMesh.lookupObject<fieldType>(fieldName_);
-        }
-    }
-    else
-    {
-        return nbrMesh.lookupObject<fieldType>(fieldName_);
-    }
-}
-
-
-template<class Type>
-const interpolation<Type>&
-mappedFixedValueFvPatchField<Type>::interpolator() const
+const mappedPatchBase& mappedFixedValueFvPatchField<Type>::mapper
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
 {
-    if (!interpolator_.valid())
+    if (!isA<mappedPatchBase>(p.patch()))
     {
-        interpolator_ = interpolation<Type>::New
+        FatalErrorIn
         (
-            interpolationScheme_,
-            sampleField()
-        );
+            "mappedFixedValueFvPatchField<Type>::mapper()"
+        )   << "\n    patch type '" << p.patch().type()
+            << "' not type '" << mappedPatchBase::typeName << "'"
+            << "\n    for patch " << p.patch().name()
+            << " of field " << iF.name()
+            << " in file " << iF.objectPath()
+            << exit(FatalError);
     }
-    return interpolator_();
+    return refCast<const mappedPatchBase>(p.patch());
 }
 
 
@@ -222,140 +129,7 @@ void mappedFixedValueFvPatchField<Type>::updateCoeffs()
         return;
     }
 
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    // Since we're inside initEvaluate/evaluate there might be processor
-    // comms underway. Change the tag we use.
-    int oldTag = UPstream::msgType();
-    UPstream::msgType() = oldTag + 1;
-
-    // Get the scheduling information from the mappedPatchBase
-    const mappedPatchBase& mpp = refCast<const mappedPatchBase>
-    (
-        mappedFixedValueFvPatchField<Type>::patch().patch()
-    );
-
-    const fvMesh& thisMesh = this->patch().boundaryMesh().mesh();
-    const fvMesh& nbrMesh = refCast<const fvMesh>(mpp.sampleMesh());
-
-    // Result of obtaining remote values
-    Field<Type> newValues;
-
-    switch (mpp.mode())
-    {
-        case mappedPatchBase::NEARESTCELL:
-        {
-            const mapDistribute& distMap = mpp.map();
-
-            if (interpolationScheme_ != interpolationCell<Type>::typeName)
-            {
-                // Send back sample points to the processor that holds the cell
-                vectorField samples(mpp.samplePoints());
-                distMap.reverseDistribute
-                (
-                    (mpp.sameRegion() ? thisMesh.nCells() : nbrMesh.nCells()),
-                    point::max,
-                    samples
-                );
-
-                const interpolation<Type>& interp = interpolator();
-
-                newValues.setSize(samples.size(), pTraits<Type>::max);
-                forAll(samples, cellI)
-                {
-                    if (samples[cellI] != point::max)
-                    {
-                        newValues[cellI] = interp.interpolate
-                        (
-                            samples[cellI],
-                            cellI
-                        );
-                    }
-                }
-            }
-            else
-            {
-                newValues = sampleField();
-            }
-
-            distMap.distribute(newValues);
-
-            break;
-        }
-        case mappedPatchBase::NEARESTPATCHFACE:
-        case mappedPatchBase::NEARESTPATCHFACEAMI:
-        {
-            const label nbrPatchID =
-                nbrMesh.boundaryMesh().findPatchID(mpp.samplePatch());
-
-            if (nbrPatchID < 0)
-            {
-                FatalErrorIn
-                (
-                    "void mappedFixedValueFvPatchField<Type>::updateCoeffs()"
-                )<< "Unable to find sample patch " << mpp.samplePatch()
-                 << " in region " << mpp.sampleRegion()
-                 << " for patch " << this->patch().name() << nl
-                 << abort(FatalError);
-            }
-
-            const fieldType& nbrField = sampleField();
-
-            newValues = nbrField.boundaryField()[nbrPatchID];
-            mpp.distribute(newValues);
-
-            break;
-        }
-        case mappedPatchBase::NEARESTFACE:
-        {
-            Field<Type> allValues(nbrMesh.nFaces(), pTraits<Type>::zero);
-
-            const fieldType& nbrField = sampleField();
-
-            forAll(nbrField.boundaryField(), patchI)
-            {
-                const fvPatchField<Type>& pf =
-                    nbrField.boundaryField()[patchI];
-                label faceStart = pf.patch().start();
-
-                forAll(pf, faceI)
-                {
-                    allValues[faceStart++] = pf[faceI];
-                }
-            }
-
-            mpp.distribute(allValues);
-            newValues.transfer(allValues);
-
-            break;
-        }
-        default:
-        {
-            FatalErrorIn
-            (
-                "mappedFixedValueFvPatchField<Type>::updateCoeffs()"
-            )<< "Unknown sampling mode: " << mpp.mode()
-             << nl << abort(FatalError);
-        }
-    }
-
-    if (setAverage_)
-    {
-        Type averagePsi =
-            gSum(this->patch().magSf()*newValues)
-           /gSum(this->patch().magSf());
-
-        if (mag(averagePsi)/mag(average_) > 0.5)
-        {
-            newValues *= mag(average_)/mag(averagePsi);
-        }
-        else
-        {
-            newValues += (average_ - averagePsi);
-        }
-    }
-
-    this->operator==(newValues);
+    this->operator==(this->mappedField());
 
     if (debug)
     {
@@ -368,9 +142,6 @@ void mappedFixedValueFvPatchField<Type>::updateCoeffs()
             << endl;
     }
 
-    // Restore tag
-    UPstream::msgType() = oldTag;
-
     fixedValueFvPatchField<Type>::updateCoeffs();
 }
 
@@ -379,11 +150,7 @@ template<class Type>
 void mappedFixedValueFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
-    os.writeKeyword("fieldName") << fieldName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl;
-    os.writeKeyword("average") << average_ << token::END_STATEMENT << nl;
-    os.writeKeyword("interpolationScheme") << interpolationScheme_
-        << token::END_STATEMENT << nl;
+    mappedPatchFieldBase<Type>::write(os);
     this->writeEntry("value", os);
 }
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.H
index dd52f2a2b635b725368edc81cf7654ccdeab7485..6267e88ab228791cf362cdeeb494351f7e5653fc 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedValue/mappedFixedValueFvPatchField.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
@@ -83,7 +83,8 @@ SourceFiles
 #define mappedFixedValueFvPatchField_H
 
 #include "fixedValueFvPatchFields.H"
-#include "interpolation.H"
+//#include "interpolation.H"
+#include "mappedPatchFieldBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -97,37 +98,19 @@ namespace Foam
 template<class Type>
 class mappedFixedValueFvPatchField
 :
-    public fixedValueFvPatchField<Type>
+    public fixedValueFvPatchField<Type>,
+    public mappedPatchFieldBase<Type>
 {
 
 protected:
 
-    // Protected data
-
-        //- Name of field to sample - defaults to field associated with this
-        //  patchField if not specified
-        word fieldName_;
-
-        //- If true adjust the mapped field to maintain average value average_
-        const bool setAverage_;
-
-        //- Average value the mapped field is adjusted to maintain if
-        //  setAverage_ is set true
-        const Type average_;
-
-        //- Interpolation scheme to use for nearestcell mode
-        word interpolationScheme_;
-
-        mutable autoPtr<interpolation<Type> > interpolator_;
-
-
     // Protected Member Functions
 
-        //- Field to sample. Either on my or nbr mesh
-        const GeometricField<Type, fvPatchField, volMesh>& sampleField() const;
-
-        //- Access the interpolation method
-        const interpolation<Type>& interpolator() const;
+        const mappedPatchBase& mapper
+        (
+            const fvPatch& p,
+            const DimensionedField<Type, volMesh>& iF
+        );
 
 
 public:
diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C
index f5622451cbf45706cc3c7569aff24b2ea9bb2d1b..fa3d6d7a315e0fcea6297366534e853e408fa9bc 100644
--- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C
@@ -198,10 +198,61 @@ void Foam::fvPatchField<Type>::patchInternalField(Field<Type>& pif) const
 template<class Type>
 void Foam::fvPatchField<Type>::autoMap
 (
-    const fvPatchFieldMapper& m
+    const fvPatchFieldMapper& mapper
 )
 {
-    Field<Type>::autoMap(m);
+    Field<Type>& f = *this;
+
+    if (!this->size())
+    {
+        f.setSize(mapper.size());
+        if (f.size())
+        {
+            f = this->patchInternalField();
+        }
+    }
+    else
+    {
+        // Map all faces provided with mapping data
+        Field<Type>::autoMap(mapper);
+
+        // For unmapped faces set to internal field value (zero-gradient)
+        if
+        (
+            mapper.direct()
+         && &mapper.directAddressing()
+         && mapper.directAddressing().size()
+        )
+        {
+            Field<Type> pif(this->patchInternalField());
+
+            const labelList& mapAddressing = mapper.directAddressing();
+
+            forAll(mapAddressing, i)
+            {
+                if (mapAddressing[i] < 0)
+                {
+                    f[i] = pif[i];
+                }
+            }
+        }
+        else if (!mapper.direct() && mapper.addressing().size())
+        {
+            Field<Type> pif(this->patchInternalField());
+
+            const labelListList& mapAddressing = mapper.addressing();
+
+            forAll(mapAddressing, i)
+            {
+                const labelList& localAddrs = mapAddressing[i];
+
+                if (!localAddrs.size())
+                {
+                    f[i] = pif[i];
+                }
+            }
+        }
+    }
 }
 
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
index d1ed7e01c0489d38abfa66be1a1e1b57f9a248b4..1a699590953e2e94da8d701e8f989521169ebf47 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -171,7 +171,10 @@ surfaceInterpolationScheme<Type>::interpolate
                "(const GeometricField<Type, fvPatchField, volMesh>&, "
                "const tmp<surfaceScalarField>&, "
                "const tmp<surfaceScalarField>&) : "
-               "interpolating volTypeField from cells to faces "
+               "interpolating "
+            << vf.type() << " "
+            << vf.name()
+            << " from cells to faces "
                "without explicit correction"
             << endl;
     }
@@ -252,7 +255,10 @@ surfaceInterpolationScheme<Type>::interpolate
         Info<< "surfaceInterpolationScheme<Type>::interpolate"
                "(const GeometricField<Type, fvPatchField, volMesh>&, "
                "const tmp<surfaceScalarField>&) : "
-               "interpolating volTypeField from cells to faces "
+               "interpolating "
+            << vf.type() << " "
+            << vf.name()
+            << " from cells to faces "
                "without explicit correction"
             << endl;
     }
@@ -326,7 +332,10 @@ surfaceInterpolationScheme<Type>::interpolate
     {
         Info<< "surfaceInterpolationScheme<Type>::interpolate"
                "(const GeometricField<Type, fvPatchField, volMesh>&) : "
-            << "interpolating volTypeField from cells to faces"
+               "interpolating "
+            << vf.type() << " "
+            << vf.name()
+            << " from cells to faces"
             << endl;
     }
 
diff --git a/src/meshTools/sets/topoSets/faceZoneSet.C b/src/meshTools/sets/topoSets/faceZoneSet.C
index 10d9c0b9bfd0887e88a67969651c880822601cba..5b4849c1eff51e974aa474be8bc2b47e75e71691 100644
--- a/src/meshTools/sets/topoSets/faceZoneSet.C
+++ b/src/meshTools/sets/topoSets/faceZoneSet.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,8 +49,8 @@ void faceZoneSet::updateSet()
 {
     labelList order;
     sortedOrder(addressing_, order);
-    inplaceReorder(order, addressing_);
-    inplaceReorder(order, flipMap_);
+    addressing_ = UIndirectList<label>(addressing_, order)();
+    flipMap_ = UIndirectList<bool>(flipMap_, order)();
 
     faceSet::clearStorage();
     faceSet::resize(2*addressing_.size());
diff --git a/src/thermophysicalModels/specie/transport/const/constTransport.H b/src/thermophysicalModels/specie/transport/const/constTransport.H
index c1e3120000c6a541fec680e8da07e744fddd1276..ec6bb6573732b1890023bfed41b23a0a06e6ba19 100644
--- a/src/thermophysicalModels/specie/transport/const/constTransport.H
+++ b/src/thermophysicalModels/specie/transport/const/constTransport.H
@@ -161,10 +161,13 @@ public:
 
     // Member operators
 
-        inline constTransport& operator=
-        (
-            const constTransport&
-        );
+        inline constTransport& operator=(const constTransport&);
+
+        inline void operator+=(const constTransport&);
+
+        inline void operator-=(const constTransport&);
+
+        inline void operator*=(const scalar);
 
 
     // Friend operators
diff --git a/src/thermophysicalModels/specie/transport/const/constTransportI.H b/src/thermophysicalModels/specie/transport/const/constTransportI.H
index e58aa79efa95828f2525b1349de8360ae48fcb0b..1581de144ef9d8324272aeb5d1fddb7985d7cff2 100644
--- a/src/thermophysicalModels/specie/transport/const/constTransportI.H
+++ b/src/thermophysicalModels/specie/transport/const/constTransportI.H
@@ -143,6 +143,52 @@ inline Foam::constTransport<Thermo>& Foam::constTransport<Thermo>::operator=
 }
 
 
+template<class Thermo>
+inline void Foam::constTransport<Thermo>::operator+=
+(
+    const constTransport<Thermo>& st
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator+=(st);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = st.nMoles()/this->nMoles();
+
+    mu_ = molr1*mu_ + molr2*st.mu_;
+    rPr_ = molr1*rPr_ + molr2*st.rPr_;
+}
+
+
+template<class Thermo>
+inline void Foam::constTransport<Thermo>::operator-=
+(
+    const constTransport<Thermo>& st
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator-=(st);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = st.nMoles()/this->nMoles();
+
+    mu_ = molr1*mu_ - molr2*st.mu_;
+    rPr_ = molr1*rPr_ - molr2*st.rPr_;
+}
+
+
+template<class Thermo>
+inline void Foam::constTransport<Thermo>::operator*=
+(
+    const scalar s
+)
+{
+    Thermo::operator*=(s);
+}
+
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 template<class Thermo>
diff --git a/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransport.H b/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransport.H
index f8ac0722dd6178efda9d62653d2c218c1b68c1fd..55f5f27e5a27bc0222b83bd2d5f2386f874e1ec0 100644
--- a/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransport.H
+++ b/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransport.H
@@ -180,10 +180,13 @@ public:
 
     // Member operators
 
-        inline sutherlandTransport& operator=
-        (
-            const sutherlandTransport&
-        );
+        inline sutherlandTransport& operator=(const sutherlandTransport&);
+
+        inline void operator+=(const sutherlandTransport&);
+
+        inline void operator-=(const sutherlandTransport&);
+
+        inline void operator*=(const scalar);
 
 
     // Friend operators
diff --git a/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H b/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H
index ac79418ed34435d4ca2c16857ef03316a5c8d897..56fcaf6e8f080530b2f4d96be0914360daf92bfe 100644
--- a/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H
+++ b/src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H
@@ -180,6 +180,52 @@ Foam::sutherlandTransport<Thermo>::operator=
 }
 
 
+template<class Thermo>
+inline void Foam::sutherlandTransport<Thermo>::operator+=
+(
+    const sutherlandTransport<Thermo>& st
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator+=(st);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = st.nMoles()/this->nMoles();
+
+    As_ = molr1*As_ + molr2*st.As_;
+    Ts_ = molr1*Ts_ + molr2*st.Ts_;
+}
+
+
+template<class Thermo>
+inline void Foam::sutherlandTransport<Thermo>::operator-=
+(
+    const sutherlandTransport<Thermo>& st
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator-=(st);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = st.nMoles()/this->nMoles();
+
+    As_ = molr1*As_ - molr2*st.As_;
+    Ts_ = molr1*Ts_ - molr2*st.Ts_;
+}
+
+
+template<class Thermo>
+inline void Foam::sutherlandTransport<Thermo>::operator*=
+(
+    const scalar s
+)
+{
+    Thermo::operator*=(s);
+}
+
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 template<class Thermo>
diff --git a/src/triSurface/triSurface/interfaces/STL/readSTLASCII.L b/src/triSurface/triSurface/interfaces/STL/readSTLASCII.L
index f2fb5c6e4d2975fc92047e28549868b5d3298a58..3c732c1a8d7088097fa6f6b2945cc46c4e411d08 100644
--- a/src/triSurface/triSurface/interfaces/STL/readSTLASCII.L
+++ b/src/triSurface/triSurface/interfaces/STL/readSTLASCII.L
@@ -55,7 +55,7 @@ int yyFlexLexer::yylex()
 // It is called by yylex but is not used as the mechanism to change file.
 // See <<EOF>>
 //! \cond dummy
-#if YY_FLEX_SUBMINOR_VERSION < 34 || YY_FLEX_SUBMINOR_VERSION > 36
+#if YY_FLEX_SUBMINOR_VERSION < 34
 extern "C" int yywrap()
 #else
 int yyFlexLexer::yywrap()
diff --git a/src/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.C b/src/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.C
index 060ed5bda3c37385932ff2f3d66c33661056d263..08ce6197c6543cd4a14c4fc32790b4475091290b 100644
--- a/src/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.C
+++ b/src/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.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
@@ -66,8 +66,12 @@ Foam::anisotropicFilter::anisotropicFilter
         coeff_.internalField().replace
         (
             d,
-            (2.0/widthCoeff_)*mesh.V()
-           /fvc::surfaceSum(mag(mesh.Sf().component(d)))().internalField()
+            (1/widthCoeff_)*
+            sqr
+            (
+                2.0*mesh.V()
+               /fvc::surfaceSum(mag(mesh.Sf().component(d)))().internalField()
+            )
         );
     }
 }
@@ -99,8 +103,12 @@ Foam::anisotropicFilter::anisotropicFilter
         coeff_.internalField().replace
         (
             d,
-            (2.0/widthCoeff_)*mesh.V()
-            /fvc::surfaceSum(mag(mesh.Sf().component(d)))().internalField()
+            (1/widthCoeff_)*
+            sqr
+            (
+                2.0*mesh.V()
+               /fvc::surfaceSum(mag(mesh.Sf().component(d)))().internalField()
+            )
         );
     }
 }
diff --git a/src/turbulenceModels/LES/LESfilters/laplaceFilter/laplaceFilter.C b/src/turbulenceModels/LES/LESfilters/laplaceFilter/laplaceFilter.C
index 231ae2737fa9023a643ef57b2f8606652ba21343..24db120ffa99b37c454e06e473b8780f30c29ffc 100644
--- a/src/turbulenceModels/LES/LESfilters/laplaceFilter/laplaceFilter.C
+++ b/src/turbulenceModels/LES/LESfilters/laplaceFilter/laplaceFilter.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,7 +57,7 @@ Foam::laplaceFilter::laplaceFilter(const fvMesh& mesh, scalar widthCoeff)
         calculatedFvPatchScalarField::typeName
     )
 {
-    coeff_.internalField() = pow(mesh.V(), 2.0/3.0)/widthCoeff_;
+    coeff_.dimensionedInternalField() = pow(mesh.V(), 2.0/3.0)/widthCoeff_;
 }
 
 
@@ -78,7 +78,7 @@ Foam::laplaceFilter::laplaceFilter(const fvMesh& mesh, const dictionary& bd)
         calculatedFvPatchScalarField::typeName
     )
 {
-    coeff_.internalField() = pow(mesh.V(), 2.0/3.0)/widthCoeff_;
+    coeff_.dimensionedInternalField() = pow(mesh.V(), 2.0/3.0)/widthCoeff_;
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
index 40c4deaf33a3d1f7b052b93fc2ba5d7844cf83ea..5e125d53108ff57fc8c63829997f440aa9bb2e1a 100644
--- a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,7 +66,7 @@ volScalarField dynOneEqEddy::ck
 
     const volSymmTensorField MM
     (
-        simpleFilter_(-2.0*delta()*pow(KK, 0.5)*filter_(D))
+        simpleFilter_(-2.0*delta()*sqrt(KK)*filter_(D))
     );
 
     const volScalarField ck
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes
index 469090b5e0c522c21d90d80958ba88389a249b6d..3582eeafa1a910f67c7fc1e7194b8aee21994754 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes
@@ -27,14 +27,17 @@ gradSchemes
 
 divSchemes
 {
-    div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression 1;
-    div(phid1,p_rgh) Gauss upwind;
-    div(phid2,p_rgh) Gauss upwind;
+
+    div(rho*phi,U)  Gauss upwind;
+    div(phi,thermo:rhowater) Gauss upwind;
+    div(phi,thermo:rhoair) Gauss upwind;
     div(rho*phi,T)  Gauss upwind;
     div(rho*phi,K)  Gauss upwind;
-    div(phi,k)      Gauss vanLeer;
+    div(phi,p)      Gauss upwind;
+    div(phi,k)      Gauss upwind;
+
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution
index 9066c1c700c4ba078eae024acd15d76d07f45f55..7577b94f2117c0747b8af59a35a1d46738d3b7ab 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution
@@ -103,9 +103,9 @@ solvers
 PIMPLE
 {
     momentumPredictor no;
-    transSonic      no;
-    nOuterCorrectors 3;
-    nCorrectors     1;
+    transonic       no;
+    nOuterCorrectors 1;
+    nCorrectors     2;
     nNonOrthogonalCorrectors 0;
     nAlphaCorr      1;
     nAlphaSubCycles 1;
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/constant/transportProperties b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/constant/transportProperties
deleted file mode 100644
index 564df56f2083c586a247b9a7276aba1b90d3a4f2..0000000000000000000000000000000000000000
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/constant/transportProperties
+++ /dev/null
@@ -1,59 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      transportProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-phases (water air);
-
-water
-{
-    transportModel  Newtonian;
-    nu              1e-06;
-    rho             1000;
-    k               0; // 0.613;
-    Cv              4179;
-
-    equationOfState
-    {
-        type            perfectFluid;
-
-        rho0            1000;
-        R               3000;
-    }
-}
-
-air
-{
-    transportModel  Newtonian;
-    nu              1.589e-05;
-    rho             1;
-    k               0; // 2.63e-2;
-    Cv              721;
-
-    equationOfState
-    {
-        type            perfectFluid;
-
-        rho0            0;
-        R               287;
-    }
-}
-
-pMin            pMin [ 1 -1 -2 0 0 0 0 ] 10000;
-
-sigma           sigma [ 1 0 -2 0 0 0 0 ] 0.07;
-
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes
index 469090b5e0c522c21d90d80958ba88389a249b6d..3582eeafa1a910f67c7fc1e7194b8aee21994754 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes
@@ -27,14 +27,17 @@ gradSchemes
 
 divSchemes
 {
-    div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression 1;
-    div(phid1,p_rgh) Gauss upwind;
-    div(phid2,p_rgh) Gauss upwind;
+
+    div(rho*phi,U)  Gauss upwind;
+    div(phi,thermo:rhowater) Gauss upwind;
+    div(phi,thermo:rhoair) Gauss upwind;
     div(rho*phi,T)  Gauss upwind;
     div(rho*phi,K)  Gauss upwind;
-    div(phi,k)      Gauss vanLeer;
+    div(phi,p)      Gauss upwind;
+    div(phi,k)      Gauss upwind;
+
     div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSolution b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSolution
index 9066c1c700c4ba078eae024acd15d76d07f45f55..7577b94f2117c0747b8af59a35a1d46738d3b7ab 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSolution
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSolution
@@ -103,9 +103,9 @@ solvers
 PIMPLE
 {
     momentumPredictor no;
-    transSonic      no;
-    nOuterCorrectors 3;
-    nCorrectors     1;
+    transonic       no;
+    nOuterCorrectors 1;
+    nCorrectors     2;
     nNonOrthogonalCorrectors 0;
     nAlphaCorr      1;
     nAlphaSubCycles 1;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/epsilon b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/epsilon
index d5ee15eb79d080246921e64270bf55a4d1541884..bb5b3e0a73f083f93109b40e0d3659f5f1466b6d 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/epsilon
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/epsilon
@@ -29,6 +29,7 @@ boundaryField
     outlet
     {
         type            inletOutlet;
+        phi             phiwater;
         inletValue      uniform 0.1;
         value           uniform 0.1;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/k b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/k
index a791cb824f67debc5017f8efa4da5b3e4349556d..9701bd625c5d07a6522baa8bdfbbc7b42ba0db3c 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/k
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/k
@@ -29,6 +29,7 @@ boundaryField
     outlet
     {
         type            inletOutlet;
+        phi             phiwater;
         inletValue      uniform 1e-8;
         value           uniform 1e-8;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
index 3335348f056c638bdbaa92ea65e7fdce7fd887d7..c0ea375341437d614a1cad16c12827f45ec79cdd 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
     "div\(phi.*,U.*\)"      Gauss limitedLinearV 1;
     "div\(\(alpha.*Rc\)\)"  Gauss linear;
-    "div\(phid.*,p\)"       Gauss upwind;
+    "div\(phi.*,.*rho.*\)"  Gauss linear;
 
     "div\(alphaPhi.*,h.*\)" Gauss limitedLinear 1;
     "div\(alphaPhi.*,K.*\)" Gauss limitedLinear 1;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/0/epsilon b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/0/epsilon
index ee84cd44265e7dcde3de56307d9e9d6f7d737883..69b3e74c838b2f693281b70c870742448faeb9b8 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/0/epsilon
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/0/epsilon
@@ -31,6 +31,7 @@ boundaryField
     outlet
     {
         type               inletOutlet;
+        phi                phi2;
         inletValue         uniform 10.0;
         value              uniform 10.0;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/0/k b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/0/k
index 45302e2a18f632b4bcc79ff4bd9535e971dbbf5f..1ae05e468e21b3ba28c4aecccfafa12962e30bce 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/0/k
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/0/k
@@ -31,6 +31,7 @@ boundaryField
     outlet
     {
         type               inletOutlet;
+        phi                phi2;
         inletValue         uniform 1.0;
         value              uniform 1.0;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes
index 7073eebd7f4840dedd47c8e7e5561392ed125b02..6cbee38f7871c73a58663948be4905135d9e649e 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     "div\(alphaPhi.,U.\)"   Gauss limitedLinearV 1;
     "div\(phi.,U.\)"        Gauss limitedLinearV 1;
     "div\(\(alpha.*Rc\)\)"  Gauss linear;
-    "div\(phid.,p\)"        Gauss upwind;
+    "div\(phi.*,.*rho.*\)"  Gauss linear;
 
     "div\(alphaPhi.,(h|e).\)"   Gauss limitedLinear 1;
     "div\(alphaPhi.,(K.|p)\)"   Gauss limitedLinear 1;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/fvSchemes
index 6af6a5df121a97ccc62e0bb663ccf1331185197b..275907488b4c204290759f85a7003bb66f553902 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     "div\(alphaPhi.,U.\)"   Gauss limitedLinearV 1;
     "div\(phi.,U.\)"        Gauss limitedLinearV 1;
     "div\(\(alpha.*Rc\)\)"  Gauss linear;
-    "div\(phid.,p\)"        Gauss linear;
+    "div\(phi.*,.*rho.*\)"  Gauss linear;
 
     "div\(alphaPhi.,h.\)"   Gauss limitedLinear 1;
     "div\(alphaPhi.,K.\)"   Gauss limitedLinear 1;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed/0/Theta b/tutorials/multiphase/twoPhaseEulerFoam/bed/0/Theta
index 3a4a3b6744614cd730f8805ad21820c9e99117e8..2ac308e7d1ce434b094e84b917904db8fd356bb7 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed/0/Theta
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed/0/Theta
@@ -30,6 +30,7 @@ boundaryField
     top
     {
         type                inletOutlet;
+        phi                 phi1;
         inletValue          uniform 1.0e-8;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed/0/epsilon b/tutorials/multiphase/twoPhaseEulerFoam/bed/0/epsilon
index 79b3e97bad66442ae9f9dda9d0eeaaa0edfa2ad3..79fe0eda3fb27eb2a2c43092f0df67f3bb7f12e4 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed/0/epsilon
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed/0/epsilon
@@ -29,6 +29,7 @@ boundaryField
     top
     {
         type               inletOutlet;
+        phi                phi2;
         inletValue         uniform 10.0;
         value              uniform 10.0;
     }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed/0/k b/tutorials/multiphase/twoPhaseEulerFoam/bed/0/k
index ebc748ced17d51773bdff65585a2e0a278dcd022..188184e8b4dbf16123ec79623d9d19789baa973f 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed/0/k
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed/0/k
@@ -29,6 +29,7 @@ boundaryField
     top
     {
         type               inletOutlet;
+        phi                phi2;
         inletValue         uniform 1.0;
         value              uniform 1.0;
     }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed2/0/epsilon b/tutorials/multiphase/twoPhaseEulerFoam/bed2/0/epsilon
index ee88570d3991e180c4cc884e4202b4c61309d80f..4ff943aac34460251e6bdbe47cc276215d1f766e 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed2/0/epsilon
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed2/0/epsilon
@@ -27,6 +27,7 @@ boundaryField
     outlet
     {
         type               inletOutlet;
+        phi                phi2;
         inletValue         uniform 10.0;
         value              uniform 10.0;
     }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed2/0/k b/tutorials/multiphase/twoPhaseEulerFoam/bed2/0/k
index 0d6e80afbab447df5016caddee63f112c9ce03b5..6b3ddf98cffd1674871f61e8f5f47d38d23fcbc8 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed2/0/k
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed2/0/k
@@ -27,6 +27,7 @@ boundaryField
     outlet
     {
         type               inletOutlet;
+        phi                phi2;
         inletValue         uniform 1.0;
         value              uniform 1.0;
     }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/Theta b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/Theta
index e45304b83464ea6a9568531b35b570abae0d768f..523be90f325771799ebe37548e2fb4e5e8b460e9 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/Theta
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/Theta
@@ -29,6 +29,7 @@ boundaryField
     outlet
     {
         type            inletOutlet;
+        phi             phi1;
         inletValue      uniform 1.0e-7;
         value           uniform 1.0e-7;
     }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/epsilon b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/epsilon
index d5ee15eb79d080246921e64270bf55a4d1541884..d54d0a402bae4eb3e2e8f7916d06fb9de239be57 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/epsilon
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/epsilon
@@ -29,6 +29,7 @@ boundaryField
     outlet
     {
         type            inletOutlet;
+        phi             phi2;
         inletValue      uniform 0.1;
         value           uniform 0.1;
     }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/k b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/k
index a791cb824f67debc5017f8efa4da5b3e4349556d..bf8d7cc467cd6a46871954606b86a3132f02972a 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/k
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/0/k
@@ -29,6 +29,7 @@ boundaryField
     outlet
     {
         type            inletOutlet;
+        phi             phi2;
         inletValue      uniform 1e-8;
         value           uniform 1e-8;
     }