diff --git a/applications/solvers/multiphase/cavitatingFoam/CourantNo.H b/applications/solvers/multiphase/cavitatingFoam/CourantNo.H
index 02a6c86db7..a5bae74528 100644
--- a/applications/solvers/multiphase/cavitatingFoam/CourantNo.H
+++ b/applications/solvers/multiphase/cavitatingFoam/CourantNo.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ if (mesh.nInternalFaces())
 {
     scalarField sumPhi
     (
-        fvc::surfaceSum(mag(phiv))().internalField()
+        fvc::surfaceSum(mag(phi))().internalField()
     );
 
     CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
@@ -54,7 +54,7 @@ if (mesh.nInternalFaces())
     )*runTime.deltaTValue();
 }
 
-Info<< "phiv Courant Number mean: " << meanCoNum
+Info<< "phi Courant Number mean: " << meanCoNum
     << " max: " << CoNum
     << " acoustic max: " << acousticCoNum
     << endl;
diff --git a/applications/solvers/multiphase/cavitatingFoam/UEqn.H b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
index 1ca0b9f0a5..8edf0d3934 100644
--- a/applications/solvers/multiphase/cavitatingFoam/UEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
@@ -1,7 +1,7 @@
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
-      + fvm::div(phi, U)
+      + fvm::div(rhoPhi, U)
       + turbulence->divDevRhoReff(rho, U)
     );
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
index c3fde6b3a1..284ca97bd4 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createUf.H"
     #include "createPcorrTypes.H"
-    #include "compressibleCourantNo.H"
+    #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -86,12 +86,12 @@ int main(int argc, char *argv[])
             if (correctPhi)
             {
                 // Calculate absolute flux from the mapped surface velocity
-                phiv = mesh.Sf() & Uf;
+                phi = mesh.Sf() & Uf;
 
                 #include "correctPhi.H"
 
                 // Make the flux relative to the mesh motion
-                fvc::makeRelative(phiv, U);
+                fvc::makeRelative(phi, U);
             }
         }
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
index c05d664838..4bb1ce9d0d 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
@@ -14,7 +14,7 @@ if (mesh.changing())
         {
             U.boundaryField()[patchI].evaluate();
 
-            phiv.boundaryField()[patchI] =
+            phi.boundaryField()[patchI] =
                 U.boundaryField()[patchI]
               & mesh.Sf().boundaryField()[patchI];
         }
@@ -37,21 +37,21 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    surfaceScalarField rhof(fvc::interpolate(rho, "div(phiv,rho)"));
+    surfaceScalarField rhof(fvc::interpolate(rho, "div(phi,rho)"));
     dimensionedScalar rAUf("rAUf", dimTime, 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(rAUf, pcorr) == fvc::ddt(rho) + fvc::div(phiv*rhof)
+            fvm::laplacian(rAUf, pcorr) == fvc::ddt(rho) + fvc::div(phi*rhof)
         );
 
         pcorrEqn.solve();
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phiv -= pcorrEqn.flux()/rhof;
+            phi -= pcorrEqn.flux()/rhof;
         }
     }
 }
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index a0e0e6c74d..b07846b6de 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -17,13 +17,13 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    phiv = (fvc::interpolate(HbyA) & mesh.Sf())
+    phi = (fvc::interpolate(HbyA) & mesh.Sf())
          + rhorAUf*fvc::ddtCorr(U, Uf);
-    fvc::makeRelative(phiv, U);
+    fvc::makeRelative(phi, U);
 
     surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
 
-    phiv -= phiGradp/rhof;
+    phi -= phiGradp/rhof;
 
     volScalarField rho0(rho - psi*p);
 
@@ -33,7 +33,7 @@
         (
             fvc::ddt(rho)
           + psi*correction(fvm::ddt(p))
-          + fvc::div(phiv, rho)
+          + fvc::div(phi, rho)
           + fvc::div(phiGradp)
           - fvm::laplacian(rhorAUf, p)
         );
@@ -42,7 +42,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phiv += (phiGradp + pEqn.flux())/rhof;
+            phi += (phiGradp + pEqn.flux())/rhof;
         }
     }
 
@@ -85,6 +85,6 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += n*(phiv/mesh.magSf() - (n & Uf));
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
     }
 }
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
index c4453df391..2b04d58fb5 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
     #include "readControls.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
-    #include "compressibleCourantNo.H"
+    #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
     pimpleControl pimple(mesh);
diff --git a/applications/solvers/multiphase/cavitatingFoam/createFields.H b/applications/solvers/multiphase/cavitatingFoam/createFields.H
index 605c0deece..2c834dc511 100644
--- a/applications/solvers/multiphase/cavitatingFoam/createFields.H
+++ b/applications/solvers/multiphase/cavitatingFoam/createFields.H
@@ -39,12 +39,25 @@
         mesh
     );
 
-    #include "createPhiv.H"
-    #include "compressibleCreatePhi.H"
+    #include "createPhi.H"
+
+    // Mass flux (corrected by rhoEqn.H)
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rhoPhi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        fvc::interpolate(rho)*phi
+    );
 
     Info<< "Reading transportProperties\n" << endl;
 
-    incompressibleTwoPhaseMixture twoPhaseProperties(U, phiv);
+    incompressibleTwoPhaseMixture twoPhaseProperties(U, phi);
 
     volScalarField& alphav(twoPhaseProperties.alpha1());
     alphav.oldTime();
@@ -72,5 +85,5 @@
     // Create incompressible turbulence model
     autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::turbulenceModel::New(U, phiv, twoPhaseProperties)
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index 210a30fd14..da8047914c 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -17,12 +17,12 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + rhorAUf*fvc::ddtCorr(U, phiv);
+    phi = (fvc::interpolate(HbyA) & mesh.Sf())
+         + rhorAUf*fvc::ddtCorr(U, phi);
 
     surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
 
-    phiv -= phiGradp/rhof;
+    phi -= phiGradp/rhof;
 
     while (pimple.correctNonOrthogonal())
     {
@@ -30,7 +30,7 @@
         (
             fvm::ddt(psi, p)
           - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(alphav) - pSat*fvc::ddt(psi)
-          + fvc::div(phiv, rho)
+          + fvc::div(phi, rho)
           + fvc::div(phiGradp)
           - fvm::laplacian(rhorAUf, p)
         );
@@ -39,7 +39,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phiv += (phiGradp + pEqn.flux())/rhof;
+            phi += (phiGradp + pEqn.flux())/rhof;
         }
     }
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/rhoEqn.H b/applications/solvers/multiphase/cavitatingFoam/rhoEqn.H
index d0bd6e1dad..a01823ad49 100644
--- a/applications/solvers/multiphase/cavitatingFoam/rhoEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/rhoEqn.H
@@ -2,12 +2,12 @@
     fvScalarMatrix rhoEqn
     (
         fvm::ddt(rho)
-      + fvm::div(phiv, rho)
+      + fvm::div(phi, rho)
     );
 
     rhoEqn.solve();
 
-    phi = rhoEqn.flux();
+    rhoPhi = rhoEqn.flux();
 
     Info<< "max-min rho: " << max(rho).value()
         << " " << min(rho).value() << endl;
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 66df09e5e9..49958bb34f 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -256,7 +256,8 @@ castellatedMeshControls
 
     // Whether any faceZones (as specified in the refinementSurfaces)
     // are only on the boundary of corresponding cellZones or also allow
-    // free-standing zone faces. Not used if there are no faceZones.
+    // free-standing zone faces and zone faces on boundaries. Not used if
+    // there are no faceZones.
     allowFreeStandingZoneFaces true;
 
 
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 666293c558..9db171db69 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anispulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,6 +46,9 @@ Description
 #include "zeroGradientFvPatchFields.H"
 #include "CuthillMcKeeRenumber.H"
 #include "fvMeshSubset.H"
+#include "cellSet.H"
+#include "faceSet.H"
+#include "pointSet.H"
 
 #ifdef FOAM_USE_ZOLTAN
 #   include "zoltanRenumber.H"
@@ -703,6 +706,7 @@ int main(int argc, char *argv[])
     bool writeMaps = false;
     bool orderPoints = false;
     label blockSize = 0;
+    bool renumberSets = true;
 
     // Construct renumberMethod
     autoPtr<IOdictionary> renumberDictPtr;
@@ -720,7 +724,6 @@ int main(int argc, char *argv[])
 
         renumberPtr = renumberMethod::New(renumberDict);
 
-
         sortCoupledFaceCells = renumberDict.lookupOrDefault
         (
             "sortCoupledFaceCells",
@@ -763,6 +766,8 @@ int main(int argc, char *argv[])
             Info<< "Writing renumber maps (new to old) to polyMesh." << nl
                 << endl;
         }
+
+        renumberSets = renumberDict.lookupOrDefault("renumberSets", true);
     }
     else
     {
@@ -864,6 +869,54 @@ int main(int argc, char *argv[])
     PtrList<surfaceTensorField> stFlds;
     ReadFields(mesh, objects, stFlds);
 
+    // Read sets
+    PtrList<cellSet> cellSets;
+    PtrList<faceSet> faceSets;
+    PtrList<pointSet> pointSets;
+    if (renumberSets)
+    {
+        // Read sets
+        IOobjectList objects(mesh, mesh.facesInstance(), "polyMesh/sets");
+        {
+            IOobjectList cSets(objects.lookupClass(cellSet::typeName));
+            if (cSets.size())
+            {
+                Info<< "Reading cellSets:" << endl;
+                forAllConstIter(IOobjectList, cSets, iter)
+                {
+                    cellSets.append(new cellSet(*iter()));
+                    Info<< "    " << cellSets.last().name() << endl;
+                }
+            }
+        }
+        {
+            IOobjectList fSets(objects.lookupClass(faceSet::typeName));
+            if (fSets.size())
+            {
+                Info<< "Reading faceSets:" << endl;
+                forAllConstIter(IOobjectList, fSets, iter)
+                {
+                    faceSets.append(new faceSet(*iter()));
+                    Info<< "    " << faceSets.last().name() << endl;
+                }
+            }
+        }
+        {
+            IOobjectList pSets(objects.lookupClass(pointSet::typeName));
+            if (pSets.size())
+            {
+                Info<< "Reading pointSets:" << endl;
+                forAllConstIter(IOobjectList, pSets, iter)
+                {
+                    pointSets.append(new pointSet(*iter()));
+                    Info<< "    " << pointSets.last().name() << endl;
+                }
+            }
+        }
+    }
+
+
+
     Info<< endl;
 
     // From renumbering:
@@ -1058,7 +1111,6 @@ int main(int argc, char *argv[])
     mesh.updateMesh(map);
 
     // Update proc maps
-    if (cellProcAddressing.headerOk())
     if
     (
         cellProcAddressing.headerOk()
@@ -1073,7 +1125,6 @@ int main(int argc, char *argv[])
             UIndirectList<label>(cellProcAddressing, map().cellMap())
         );
     }
-    if (faceProcAddressing.headerOk())
     if
     (
         faceProcAddressing.headerOk()
@@ -1104,7 +1155,6 @@ int main(int argc, char *argv[])
             }
         }
     }
-    if (pointProcAddressing.headerOk())
     if
     (
         pointProcAddressing.headerOk()
@@ -1228,7 +1278,6 @@ int main(int argc, char *argv[])
     Info<< "Writing mesh to " << mesh.facesInstance() << endl;
 
     mesh.write();
-    if (cellProcAddressing.headerOk())
     if
     (
         cellProcAddressing.headerOk()
@@ -1238,7 +1287,6 @@ int main(int argc, char *argv[])
         cellProcAddressing.instance() = mesh.facesInstance();
         cellProcAddressing.write();
     }
-    if (faceProcAddressing.headerOk())
     if
     (
         faceProcAddressing.headerOk()
@@ -1248,7 +1296,6 @@ int main(int argc, char *argv[])
         faceProcAddressing.instance() = mesh.facesInstance();
         faceProcAddressing.write();
     }
-    if (pointProcAddressing.headerOk())
     if
     (
         pointProcAddressing.headerOk()
@@ -1258,7 +1305,6 @@ int main(int argc, char *argv[])
         pointProcAddressing.instance() = mesh.facesInstance();
         pointProcAddressing.write();
     }
-    if (boundaryProcAddressing.headerOk())
     if
     (
         boundaryProcAddressing.headerOk()
@@ -1269,7 +1315,6 @@ int main(int argc, char *argv[])
         boundaryProcAddressing.write();
     }
 
-
     if (writeMaps)
     {
         // For debugging: write out region
@@ -1334,6 +1379,28 @@ int main(int argc, char *argv[])
         ).write();
     }
 
+    if (renumberSets)
+    {
+        forAll(cellSets, i)
+        {
+            cellSets[i].updateMesh(map());
+            cellSets[i].instance() = mesh.facesInstance();
+            cellSets[i].write();
+        }
+        forAll(faceSets, i)
+        {
+            faceSets[i].updateMesh(map());
+            faceSets[i].instance() = mesh.facesInstance();
+            faceSets[i].write();
+        }
+        forAll(pointSets, i)
+        {
+            pointSets[i].updateMesh(map());
+            pointSets[i].instance() = mesh.facesInstance();
+            pointSets[i].write();
+        }
+    }
+
     Info<< "\nEnd.\n" << endl;
 
     return 0;
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
index 2cd1a68ca5..8b44fd1c4a 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
@@ -34,6 +34,8 @@ sortCoupledFaceCells false;
 // Optional entry: sort points into internal and boundary points
 //orderPoints false;
 
+// Optional: suppress renumbering cellSets,faceSets,pointSets
+//renumberSets false;
 
 
 method          CuthillMcKee;
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
index 267b8dc844..64ce70b679 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -195,41 +195,44 @@ void Foam::ensightMesh::correct()
     {
         forAll(mesh_.boundary(), patchi)
         {
-            const polyPatch& p = mesh_.boundaryMesh()[patchi];
-
-            labelList& tris = boundaryFaceSets_[patchi].tris;
-            labelList& quads = boundaryFaceSets_[patchi].quads;
-            labelList& polys = boundaryFaceSets_[patchi].polys;
-
-            tris.setSize(p.size());
-            quads.setSize(p.size());
-            polys.setSize(p.size());
-
-            label nTris = 0;
-            label nQuads = 0;
-            label nPolys = 0;
-
-            forAll(p, faceI)
+            if (mesh_.boundary()[patchi].size())
             {
-                const face& f = p[faceI];
+                const polyPatch& p = mesh_.boundaryMesh()[patchi];
 
-                if (f.size() == 3)
+                labelList& tris = boundaryFaceSets_[patchi].tris;
+                labelList& quads = boundaryFaceSets_[patchi].quads;
+                labelList& polys = boundaryFaceSets_[patchi].polys;
+
+                tris.setSize(p.size());
+                quads.setSize(p.size());
+                polys.setSize(p.size());
+
+                label nTris = 0;
+                label nQuads = 0;
+                label nPolys = 0;
+
+                forAll(p, faceI)
                 {
-                    tris[nTris++] = faceI;
-                }
-                else if (f.size() == 4)
-                {
-                    quads[nQuads++] = faceI;
-                }
-                else
-                {
-                    polys[nPolys++] = faceI;
+                    const face& f = p[faceI];
+
+                    if (f.size() == 3)
+                    {
+                        tris[nTris++] = faceI;
+                    }
+                    else if (f.size() == 4)
+                    {
+                        quads[nQuads++] = faceI;
+                    }
+                    else
+                    {
+                        polys[nPolys++] = faceI;
+                    }
                 }
+
+                tris.setSize(nTris);
+                quads.setSize(nQuads);
+                polys.setSize(nPolys);
             }
-
-            tris.setSize(nTris);
-            quads.setSize(nQuads);
-            polys.setSize(nPolys);
         }
     }
 
@@ -240,9 +243,12 @@ void Foam::ensightMesh::correct()
 
         if (patchNames_.empty() || patchNames_.found(patchName))
         {
-            nfp.nTris   = boundaryFaceSets_[patchi].tris.size();
-            nfp.nQuads  = boundaryFaceSets_[patchi].quads.size();
-            nfp.nPolys  = boundaryFaceSets_[patchi].polys.size();
+            if (mesh_.boundary()[patchi].size())
+            {
+                nfp.nTris   = boundaryFaceSets_[patchi].tris.size();
+                nfp.nQuads  = boundaryFaceSets_[patchi].quads.size();
+                nfp.nPolys  = boundaryFaceSets_[patchi].polys.size();
+            }
         }
 
         reduce(nfp.nTris, sumOp<label>());
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C b/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
index 6015f1167d..cba5433426 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -227,7 +227,8 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
                     newBoundary[patchI].start(),
                     patchI,
                     boundary_,
-                    newBoundary[patchI].type()
+                    newBoundary[patchI].physicalType(),
+                    newBoundary[patchI].inGroups()
                 );
             }
         }
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
index 643e06a479..a88b194706 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,6 +107,30 @@ Foam::polyPatch::polyPatch
 }
 
 
+Foam::polyPatch::polyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const polyBoundaryMesh& bm,
+    const word& physicalType,
+    const wordList& inGroups
+)
+:
+    patchIdentifier(name, index, physicalType, inGroups),
+    primitivePatch
+    (
+        faceSubList(bm.mesh().faces(), size, start),
+        bm.mesh().points()
+    ),
+    start_(start),
+    boundaryMesh_(bm),
+    faceCellsPtr_(NULL),
+    mePtr_(NULL)
+{}
+
+
 Foam::polyPatch::polyPatch
 (
     const word& name,
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
index 94c743f740..cb2ba4abd9 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -176,6 +176,18 @@ public:
             const word& patchType
         );
 
+        //- Construct from components
+        polyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const polyBoundaryMesh& bm,
+            const word& physicalType,
+            const wordList& inGroups
+        );
+
         //- Construct from dictionary
         polyPatch
         (
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.C
index 94084ad94f..ddc6099121 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.C
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,6 @@ addToRunTimeSelectionTable(fvPatch, wedgeFvPatch, polyPatch);
 
 // * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
 
-//- Construct from polyPatch
 wedgeFvPatch::wedgeFvPatch(const polyPatch& patch, const fvBoundaryMesh& bm)
 :
     fvPatch(patch, bm),
@@ -47,6 +46,15 @@ wedgeFvPatch::wedgeFvPatch(const polyPatch& patch, const fvBoundaryMesh& bm)
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::vectorField> Foam::wedgeFvPatch::delta() const
+{
+    const vectorField nHat(nf());
+    return nHat*(nHat & (Cf() - Cn()));
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.H
index 81f6a5f320..16fab3e42b 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.H
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,6 +83,9 @@ public:
             {
                 return wedgePolyPatch_.cellT();
             }
+
+            //- Return cell-centre to face normal vector
+            virtual tmp<vectorField> delta() const;
 };
 
 
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/0/p b/tutorials/multiphase/cavitatingFoam/les/throttle/0/p
index f484f1a42e..22afb75d57 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/0/p
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/0/p
@@ -24,7 +24,6 @@ boundaryField
     {
         type            totalPressure;
         U               U;
-        phi             phiv;
         rho             rho;
         psi             none;
         gamma           1;
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
index 6b7b46498b..17346054f6 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
@@ -28,9 +28,9 @@ interpolationSchemes
 divSchemes
 {
     default         none;
-    div(phiv,rho)   Gauss vanLeer;
-    div(phi,U)      Gauss LUST grad(U);
-    div(phiv,k)     Gauss LUST grad(k);
+    div(phi,rho)    Gauss vanLeer;
+    div(rhoPhi,U)   Gauss LUST grad(U);
+    div(phi,k)      Gauss LUST grad(k);
     div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution
index 1785af6f4e..ca2022df7b 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution
@@ -78,5 +78,13 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
+relaxationFactors
+{
+    equations
+    {
+        ".*" 1;
+    }
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/0.org/p b/tutorials/multiphase/cavitatingFoam/les/throttle3D/0.org/p
index 58cc484ff7..a85f4b9d0d 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/0.org/p
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/0.org/p
@@ -24,7 +24,6 @@ boundaryField
     {
         type            totalPressure;
         U               U;
-        phi             phiv;
         rho             rho;
         psi             none;
         gamma           1;
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
index f7e2958bf6..64aa6e54b9 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
@@ -29,9 +29,9 @@ divSchemes
 {
     default         none;
 
-    div(phiv,rho)   Gauss vanLeer;
-    div(phi,U)      Gauss LUST grad(U);
-    div(phiv,k)     Gauss LUST grad(k);
+    div(phi,rho)    Gauss vanLeer;
+    div(rhoPhi,U)   Gauss LUST grad(U);
+    div(phi,k)      Gauss LUST grad(k);
     div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution
index 7133194be0..7e1a541e0a 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution
@@ -75,5 +75,13 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
+relaxationFactors
+{
+    equations
+    {
+        ".*" 1;
+    }
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/0/p b/tutorials/multiphase/cavitatingFoam/ras/throttle/0/p
index f484f1a42e..22afb75d57 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/0/p
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/0/p
@@ -24,7 +24,6 @@ boundaryField
     {
         type            totalPressure;
         U               U;
-        phi             phiv;
         rho             rho;
         psi             none;
         gamma           1;
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
index 6f20de56af..d253239d92 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
@@ -28,10 +28,10 @@ interpolationSchemes
 divSchemes
 {
     default         none;
-    div(phiv,rho)   Gauss limitedLinear 1;
-    div(phi,U)      Gauss limitedLinearV 1;
-    div(phiv,omega) Gauss limitedLinear 1;
-    div(phiv,k)     Gauss limitedLinear 1;
+    div(phi,rho)    Gauss limitedLinear 1;
+    div(rhoPhi,U)   Gauss limitedLinearV 1;
+    div(phi,omega)  Gauss limitedLinear 1;
+    div(phi,k)      Gauss limitedLinear 1;
     div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution
index 1785af6f4e..d2ebcca130 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution
@@ -33,7 +33,6 @@ solvers
         relTol          0;
     }
 
-
     p
     {
         solver          GAMG;
@@ -78,5 +77,13 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
+relaxationFactors
+{
+    equations
+    {
+        ".*" 1;
+    }
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution
index bf1202d6ea..ca3c8a6104 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution
+++ b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution
@@ -74,9 +74,6 @@ PIMPLE
 
 relaxationFactors
 {
-    fields
-    {
-    }
     equations
     {
         ".*" 1;