diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index b097abf7451d8c2dc6c500e6dd8d91ad4d281e2e..ea3f6bdb33827e223b782909a19e965c666f4d36 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -34,6 +34,7 @@
         fvc::interpolate((1.0/rho1)*rAU1*phase1.turbulence().pPrime())
        *fvc::snGrad(alpha1)*mesh.magSf()
     );
+    phiP1.boundaryField() == 0;
 
     // Phase-2 pressure flux (e.g. due to particle-particle pressure)
     surfaceScalarField phiP2
@@ -42,6 +43,7 @@
         fvc::interpolate((1.0/rho2)*rAU2*phase2.turbulence().pPrime())
        *fvc::snGrad(alpha2)*mesh.magSf()
     );
+    phiP2.boundaryField() == 0;
 
     surfaceScalarField phiHbyA1
     (
@@ -126,9 +128,11 @@
         );
 
         pEqnComp1 =
-            fvc::ddt(rho1)
-          + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1)
-          + correction
+            (
+                fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1)
+              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+            )/rho1
+          + (alpha1/rho1)*correction
             (
                 psi1*fvm::ddt(p)
               + fvm::div(phid1, p) - fvm::Sp(fvc::div(phid1), p)
@@ -137,9 +141,11 @@
         pEqnComp1().relax();
 
         pEqnComp2 =
-            fvc::ddt(rho2)
-          + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2)
-          + correction
+            (
+                fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2)
+              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+            )/rho2
+          + (alpha2/rho2)*correction
             (
                 psi2*fvm::ddt(p)
               + fvm::div(phid2, p) - fvm::Sp(fvc::div(phid2), p)
@@ -150,12 +156,18 @@
     else
     {
         pEqnComp1 =
-            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p))
-          + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1);
+            (
+                fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1)
+              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+            )/rho1
+          + (alpha1*psi1/rho1)*correction(fvm::ddt(p));
 
         pEqnComp2 =
-            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p))
-          + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2);
+            (
+                fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2)
+              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+            )/rho2
+          + (alpha2*psi2/rho2)*correction(fvm::ddt(p));
     }
 
     // Cache p prior to solve for density update
@@ -171,11 +183,7 @@
 
         solve
         (
-            (
-                (alpha1/rho1)*pEqnComp1()
-              + (alpha2/rho2)*pEqnComp2()
-            )
-          + pEqnIncomp,
+            pEqnComp1() + pEqnComp2() + pEqnIncomp,
             mesh.solver(p.select(pimple.finalInnerIter()))
         );
 
@@ -201,8 +209,8 @@
 
             fluid.dgdt() =
             (
-                pos(alpha2)*(pEqnComp2 & p)/rho2
-              - pos(alpha1)*(pEqnComp1 & p)/rho1
+                pos(alpha2)*(pEqnComp2 & p)/max(alpha2, scalar(1e-3))
+              - pos(alpha1)*(pEqnComp1 & p)/max(alpha1, scalar(1e-3))
             );
 
             p.relax();
diff --git a/applications/test/HashingSpeed/Test-HashingSpeed.C b/applications/test/HashingSpeed/Test-HashingSpeed.C
index 1b42b5d34da58037fbcc52674960f5b385ffe9a0..c3cbda702d8d172b4884b96dfa6f9be4a5f7a88f 100644
--- a/applications/test/HashingSpeed/Test-HashingSpeed.C
+++ b/applications/test/HashingSpeed/Test-HashingSpeed.C
@@ -1016,7 +1016,6 @@ int i;
 // found somewhere in the 2nd addition
 uint32_t stroustrup (const char * s, int len) {
     uint32_t h;
-    int i;
 
     for (register int i=0; i < len; ++s, ++i)
     {
diff --git a/applications/test/List/Test-List.C b/applications/test/List/Test-List.C
index ca303d53b72f6cb3af70b19d2c6a6e172979c62e..719a5ebc5c7e7e05ceaf9e274156e31b5f05f9db 100644
--- a/applications/test/List/Test-List.C
+++ b/applications/test/List/Test-List.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
@@ -137,7 +137,7 @@ int main(int argc, char *argv[])
     Info<< "    std::list constructed from Foam::List: ";
 
     std::list<vector>::iterator it;
-    for (it=stlList.begin(); it != stlList.end(); it++)
+    for (it=stlList.begin(); it != stlList.end(); ++it)
     {
         Info<< *it << " ";
     }
diff --git a/applications/test/string/Test-string.C b/applications/test/string/Test-string.C
index 8494e8a104585e972b742a042e0f5a37761e06be..f39be8c20629f8072e1d48fd8848263b21f9f59d 100644
--- a/applications/test/string/Test-string.C
+++ b/applications/test/string/Test-string.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
@@ -79,7 +79,7 @@ int main(int argc, char *argv[])
     while (iter != iter2)
     {
         Info<< *iter;
-        iter++;
+        ++iter;
     }
     Info<< "<\n";
 
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index bf8eba6e0d6f02afd597fa73ca4f25cae914f430..9db10ea85eb8673efaeec2d12a77a1a46cf1105c 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -474,10 +474,12 @@ meshQualityControls
 // Advanced
 
 // Flags for optional output
-// 0 : only write final meshes
-// 1 : write intermediate meshes
-// 2 : write volScalarField with cellLevel for postprocessing
-// 4 : write current intersections as .obj files
+//  0 : only write final meshes
+//  1 : write intermediate meshes
+//  2 : write volScalarField with cellLevel for postprocessing
+//  4 : write current mesh intersections as .obj files
+//  8 : write information about explicit feature edge refinement
+// 16 : write information about layers
 debug 0;
 
 // Merge tolerance. Is fraction of overall bounding box of initial mesh.
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/Allrun b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/Allrun
deleted file mode 100755
index 11f1cb983e2ab9bcc09f950117ff12b2a5952e32..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/Allrun
+++ /dev/null
@@ -1,22 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
-
-# Get application name
-application=`getApplication`
-
-runApplication blockMesh
-# Make random cell numbering to get faceZone orientation random
-runApplication renumberMesh -dict system/renumberMeshDict
-
-# Construct faceZone
-runApplication topoSet
-
-runApplication decomposePar -force -cellDist
-runParallel orientFaceZone 2 f0 '(10 0 0)'
-runApplication reconstructParMesh -latestTime -mergeTol 1e-6
-
-
-# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/README.txt b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/README.txt
deleted file mode 100644
index df428773022fcae8b4139b76430acb31f74628ac..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/README.txt
+++ /dev/null
@@ -1 +0,0 @@
-2013-09-16 To test orientFaceZone. Make random orientation. Orient afterwards.
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/constant/polyMesh/blockMeshDict b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/constant/polyMesh/blockMeshDict
deleted file mode 100644
index 165a600c7b4a929aaa62ad1750bfbb9d279c2e65..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/constant/polyMesh/blockMeshDict
+++ /dev/null
@@ -1,75 +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;
-    object      blockMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-convertToMeters 0.1;
-
-vertices
-(
-    (0 0 0)
-    (1 0 0)
-    (1 1 0)
-    (0 1 0)
-    (0 0 0.1)
-    (1 0 0.1)
-    (1 1 0.1)
-    (0 1 0.1)
-);
-
-blocks
-(
-    hex (0 1 2 3 4 5 6 7) (20 20 1) simpleGrading (1 1 1)
-);
-
-edges
-(
-);
-
-boundary
-(
-    movingWall
-    {
-        type wall;
-        faces
-        (
-            (3 7 6 2)
-        );
-    }
-    fixedWalls
-    {
-        type wall;
-        faces
-        (
-            (0 4 7 3)
-            (2 6 5 1)
-            (1 5 4 0)
-        );
-    }
-    frontAndBack
-    {
-        type empty;
-        faces
-        (
-            (0 3 2 1)
-            (4 5 6 7)
-        );
-    }
-);
-
-mergePatchPairs
-(
-);
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/constant/polyMesh/boundary b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/constant/polyMesh/boundary
deleted file mode 100644
index 3d253979063ae7c46c4636ec15b41fd7ec5e7e99..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/constant/polyMesh/boundary
+++ /dev/null
@@ -1,41 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  2.2.x                                 |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       polyBoundaryMesh;
-    location    "constant/polyMesh";
-    object      boundary;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-3
-(
-    movingWall
-    {
-        type            wall;
-        nFaces          20;
-        startFace       760;
-    }
-    fixedWalls
-    {
-        type            wall;
-        nFaces          60;
-        startFace       780;
-    }
-    frontAndBack
-    {
-        type            empty;
-        inGroups        1(empty);
-        nFaces          800;
-        startFace       840;
-    }
-)
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/constant/transportProperties b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/constant/transportProperties
deleted file mode 100644
index fa1c1ca0b14b50fa06f4c8577b4998addd44b1f3..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/constant/transportProperties
+++ /dev/null
@@ -1,21 +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;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-nu              nu [ 0 2 -1 0 0 0 0 ] 0.01;
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/controlDict b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/controlDict
deleted file mode 100644
index f28aaab14691c117baa5fb158619b68487a19002..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/controlDict
+++ /dev/null
@@ -1,49 +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    "system";
-    object      controlDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-application     icoFoam;
-
-startFrom       latestTime;
-
-startTime       0;
-
-stopAt          endTime;
-
-endTime         0.5;
-
-deltaT          0.005;
-
-writeControl    timeStep;
-
-writeInterval   20;
-
-purgeWrite      0;
-
-writeFormat     ascii;
-
-writePrecision  6;
-
-writeCompression off;
-
-timeFormat      general;
-
-timePrecision   6;
-
-runTimeModifiable true;
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/decomposeParDict b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/decomposeParDict
deleted file mode 100644
index 60454a0a87e62eb6e58b3d1f1b1b06353c0369b8..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/decomposeParDict
+++ /dev/null
@@ -1,135 +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;
-    note        "mesh decomposition control dictionary";
-    object      decomposeParDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-numberOfSubdomains  2;
-
-//- Keep owner and neighbour on same processor for faces in zones:
-// preserveFaceZones (heater solid1 solid3);
-
-//- Keep owner and neighbour on same processor for faces in patches:
-//  (makes sense only for cyclic patches)
-//preservePatches (cyclic_half0 cyclic_half1);
-
-//- Keep all of faceSet on a single processor. This puts all cells
-//  connected with a point, edge or face on the same processor.
-//  (just having face connected cells might not guarantee a balanced
-//  decomposition)
-// The processor can be -1 (the decompositionMethod chooses the processor
-// for a good load balance) or explicitly provided (upsets balance).
-//singleProcessorFaceSets ((f0 -1));
-
-
-//- Use the volScalarField named here as a weight for each cell in the
-//  decomposition.  For example, use a particle population field to decompose
-//  for a balanced number of particles in a lagrangian simulation.
-// weightField dsmcRhoNMean;
-
-//method          scotch;
-method          hierarchical;
-// method          simple;
-// method          metis;
-// method          manual;
-// method          multiLevel;
-// method          structured;  // does 2D decomposition of structured mesh
-
-multiLevelCoeffs
-{
-    // Decomposition methods to apply in turn. This is like hierarchical but
-    // fully general - every method can be used at every level.
-
-    level0
-    {
-        numberOfSubdomains  64;
-        //method simple;
-        //simpleCoeffs
-        //{
-        //    n           (2 1 1);
-        //    delta       0.001;
-        //}
-        method scotch;
-    }
-    level1
-    {
-        numberOfSubdomains  4;
-        method scotch;
-    }
-}
-
-// Desired output
-
-simpleCoeffs
-{
-    n           (2 1 1);
-    delta       0.001;
-}
-
-hierarchicalCoeffs
-{
-    n           (2 1 1);
-    delta       0.001;
-    order       xyz;
-}
-
-metisCoeffs
-{
- /*
-    processorWeights
-    (
-        1
-        1
-        1
-        1
-    );
-  */
-}
-
-scotchCoeffs
-{
-    //processorWeights
-    //(
-    //    1
-    //    1
-    //    1
-    //    1
-    //);
-    //writeGraph  true;
-    //strategy "b";
-}
-
-manualCoeffs
-{
-    dataFile    "decompositionData";
-}
-
-structuredCoeffs
-{
-    // Patches to do 2D decomposition on. Structured mesh only; cells have
-    // to be in 'columns' on top of patches.
-    patches     (bottomPatch);
-}
-
-//// Is the case distributed? Note: command-line argument -roots takes
-//// precedence
-//distributed     yes;
-//// Per slave (so nProcs-1 entries) the directory above the case.
-//roots
-//(
-//    "/tmp"
-//    "/tmp"
-//);
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/fvSchemes b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/fvSchemes
deleted file mode 100644
index c311eb8961a0e5567560c8f89a87efadf31973fd..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/fvSchemes
+++ /dev/null
@@ -1,60 +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    "system";
-    object      fvSchemes;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-ddtSchemes
-{
-    default         Euler;
-}
-
-gradSchemes
-{
-    default         Gauss linear;
-    grad(p)         Gauss linear;
-}
-
-divSchemes
-{
-    default         none;
-    div(phi,U)      Gauss linear;
-}
-
-laplacianSchemes
-{
-    default         none;
-    laplacian(nu,U) Gauss linear orthogonal;
-    laplacian((1|A(U)),p) Gauss linear orthogonal;
-}
-
-interpolationSchemes
-{
-    default         linear;
-    interpolate(HbyA) linear;
-}
-
-snGradSchemes
-{
-    default         orthogonal;
-}
-
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/fvSolution b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/fvSolution
deleted file mode 100644
index cc4750f16c98b1528266fdd6ea34f2d80b412fd9..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/fvSolution
+++ /dev/null
@@ -1,46 +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    "system";
-    object      fvSolution;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-solvers
-{
-    p
-    {
-        solver          PCG;
-        preconditioner  DIC;
-        tolerance       1e-06;
-        relTol          0;
-    }
-
-    U
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-}
-
-PISO
-{
-    nCorrectors     2;
-    nNonOrthogonalCorrectors 0;
-    pRefCell        0;
-    pRefValue       0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/renumberMeshDict b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/renumberMeshDict
deleted file mode 100644
index addb04236e699fa0eb5f375e72a947688de8e634..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/renumberMeshDict
+++ /dev/null
@@ -1,110 +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;
-    note        "mesh renumbering dictionary";
-    object      renumberMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Write maps from renumbered back to original mesh
-writeMaps true;
-
-// Optional entry: sort cells on coupled boundaries to last for use with
-// e.g. nonBlockingGaussSeidel.
-sortCoupledFaceCells false;
-
-// Optional entry: renumber on a block-by-block basis. It uses a
-// blockCoeffs dictionary to construct a decompositionMethod to do
-// a block subdivision) and then applies the renumberMethod to each
-// block in turn. This can be used in large cases to keep the blocks
-// fitting in cache with all the the cache misses bunched at the end.
-// This number is the approximate size of the blocks - this gets converted
-// to a number of blocks that is the input to the decomposition method.
-//blockSize 1000;
-
-// Optional entry: sort points into internal and boundary points
-//orderPoints false;
-
-
-
-//method          CuthillMcKee;
-//method          Sloan;
-//method          manual;
-method          random;
-//method          structured;
-//method          spring;
-//method          zoltan;             // only if compiled with zoltan support
-
-//CuthillMcKeeCoeffs
-//{
-//    // Reverse CuthillMcKee (RCM) or plain
-//    reverse true;
-//}
-
-manualCoeffs
-{
-    // In system directory: new-to-original (i.e. order) labelIOList
-    dataFile "cellMap";
-}
-
-
-// For extruded (i.e. structured in one direction) meshes
-structuredCoeffs
-{
-    // Patches that mesh was extruded from. These determine the starting
-    // layer of cells
-    patches (movingWall);
-    // Method to renumber the starting layer of cells
-    method  random;
-
-    // Renumber in columns (depthFirst) or in layers
-    depthFirst true;
-
-    // Optional: reverse ordering
-    //reverse false;
-}
-
-
-springCoeffs
-{
-    // Maximum jump of cell indices. Is fraction of number of cells
-    maxCo 0.01;
-
-    // Limit the amount of movement; the fraction maxCo gets decreased
-    // with every iteration
-    freezeFraction 0.999;
-
-    // Maximum number of iterations
-    maxIter 1000;
-}
-
-
-blockCoeffs
-{
-    method          scotch;
-    //method hierarchical;
-    //hierarchicalCoeffs
-    //{
-    //    n           (1 2 1);
-    //    delta       0.001;
-    //    order       xyz;
-    //}
-}
-
-
-zoltanCoeffs
-{
-    ORDER_METHOD    LOCAL_HSFC;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/topoSetDict b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/topoSetDict
deleted file mode 100644
index fca42b45f869a6a72d3437a1ab8f7c0f9dbce513..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/topoSetDict
+++ /dev/null
@@ -1,64 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  2.2.1                                 |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      topoSetDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-actions
-(
-    // Get 'T' of faces
-    {
-        name    f0;
-        type    faceSet;
-        action  new;
-        source  boxToFace;
-        sourceInfo
-        {
-            boxes
-            (
-                (0.0499 -100 -100)(0.0501 100 100)
-                (0.03 0.0499 -100)(0.05 0.0501 100)
-            );
-        }
-    }
-
-//    // Pick up some cells
-//    {
-//        name    c0;
-//        type    cellSet;
-//        action  new;
-//        source  boxToCell;
-//        sourceInfo
-//        {
-//            boxes
-//            (
-//                (-100 -100 -100)(0.05 0.05 100)
-//                (0.05 0.05 -100)(100 100 100)
-//            );
-//        }
-//    }
-
-    // Convert faceSet to faceZone
-    {
-        name    f0;
-        type    faceZoneSet;
-        action  new;
-        source  setToFaceZone;
-        sourceInfo
-        {
-            faceSet f0;
-        }
-    }
-);
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
index 9dee4d1df94aeebde25c2bdbffd8b6e9815131d9..5f0e5c0b39df3b00054a35303fab1e9bad40f819 100644
--- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict
+++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
@@ -368,6 +368,8 @@ FoamFile
 //    {
 //        faceSet f0;       // name of faceSet
 //        cellSet c0;       // name of cellSet of slave side
+//        flip    false;    // optional: flip the faceZone (so now the cellSet
+//                          //           is the master side)
 //    }
 //
 //    // Select based on surface. Orientation from normals on surface
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.H b/applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.H
index 1f66fa276e8e0e2e5340a35c64e56fe44b029c93..54887f662a2b560fdf8417bd8b8f1b33ac88cf55 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.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
@@ -78,6 +78,26 @@ public:
             hexes(nCells),
             polys(nCells)
         {}
+
+
+    // Member Functions
+
+        void setSize(const label nCells)
+        {
+            nTets = 0;
+            nPyrs = 0;
+            nPrisms = 0;
+            nHexesWedges = 0;
+            nPolys = 0;
+
+            tets.setSize(nCells);
+            pyrs.setSize(nCells);
+            prisms.setSize(nCells);
+            wedges.setSize(nCells);
+            hexes.setSize(nCells);
+            polys.setSize(nCells);
+        }
+
 };
 
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
index 0332572eeb0d16cdb907354ac36fd9a139063803..267b8dc84421a58f13afcfa678d86befe17ffcd9 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
@@ -47,7 +47,8 @@ License
 void Foam::ensightMesh::correct()
 {
     patchPartOffset_ = 2;
-    meshCellSets_ = mesh_.nCells();
+    meshCellSets_.setSize(mesh_.nCells());
+
     boundaryFaceSets_.setSize(mesh_.boundary().size());
     allPatchNames_.clear();
     patchNames_.clear();
@@ -194,44 +195,41 @@ void Foam::ensightMesh::correct()
     {
         forAll(mesh_.boundary(), patchi)
         {
-            if (mesh_.boundary()[patchi].size())
-            {
-                const polyPatch& p = mesh_.boundaryMesh()[patchi];
+            const polyPatch& p = mesh_.boundaryMesh()[patchi];
 
-                labelList& tris = boundaryFaceSets_[patchi].tris;
-                labelList& quads = boundaryFaceSets_[patchi].quads;
-                labelList& polys = boundaryFaceSets_[patchi].polys;
+            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());
+            tris.setSize(p.size());
+            quads.setSize(p.size());
+            polys.setSize(p.size());
 
-                label nTris = 0;
-                label nQuads = 0;
-                label nPolys = 0;
+            label nTris = 0;
+            label nQuads = 0;
+            label nPolys = 0;
 
-                forAll(p, faceI)
-                {
-                    const face& f = p[faceI];
+            forAll(p, 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;
-                    }
+                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);
         }
     }
 
@@ -242,12 +240,9 @@ void Foam::ensightMesh::correct()
 
         if (patchNames_.empty() || patchNames_.found(patchName))
         {
-            if (mesh_.boundary()[patchi].size())
-            {
-                nfp.nTris   = boundaryFaceSets_[patchi].tris.size();
-                nfp.nQuads  = boundaryFaceSets_[patchi].quads.size();
-                nfp.nPolys  = boundaryFaceSets_[patchi].polys.size();
-            }
+            nfp.nTris   = boundaryFaceSets_[patchi].tris.size();
+            nfp.nQuads  = boundaryFaceSets_[patchi].quads.size();
+            nfp.nPolys  = boundaryFaceSets_[patchi].polys.size();
         }
 
         reduce(nfp.nTris, sumOp<label>());
@@ -1148,6 +1143,7 @@ void Foam::ensightMesh::write
             if (nfp.nTris || nfp.nQuads || nfp.nPolys)
             {
                 const polyPatch& p = mesh_.boundaryMesh()[patchi];
+
                 const labelList& tris = boundaryFaceSets_[patchi].tris;
                 const labelList& quads = boundaryFaceSets_[patchi].quads;
                 const labelList& polys = boundaryFaceSets_[patchi].polys;
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/PV3FoamReader/pqPV3FoamReaderPanel.cxx b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/PV3FoamReader/pqPV3FoamReaderPanel.cxx
index f795ed8eec77b749c1c94b1dd94eded37dc965dd..3b111ff59f079dc42e325a9ff565865c0438246a 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/PV3FoamReader/pqPV3FoamReaderPanel.cxx
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/PV3FoamReader/pqPV3FoamReaderPanel.cxx
@@ -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
@@ -441,8 +441,6 @@ void pqPV3FoamReaderPanel::IncludeZonesToggled()
 
 void pqPV3FoamReaderPanel::ExtrapolatePatchesToggled()
 {
-    vtkSMProperty* prop;
-
     vtkSMIntVectorProperty::SafeDownCast
     (
         this->proxy()->GetProperty("UiExtrapolatePatches")
@@ -454,8 +452,6 @@ void pqPV3FoamReaderPanel::ExtrapolatePatchesToggled()
 
 void pqPV3FoamReaderPanel::InterpolateVolFieldsToggled()
 {
-    vtkSMProperty* prop;
-
     vtkSMIntVectorProperty::SafeDownCast
     (
         this->proxy()->GetProperty("UiInterpolateVolFields")
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
index b33a9bc0918b83b0f457e4d9670da8bfbafed1ad..ff38ae0917fcbff04785468f71bfff36fae54786 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
+++ b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
@@ -358,6 +358,59 @@ Type sum(const UList<Type>& f)
 
 TMP_UNARY_FUNCTION(Type, sum)
 
+template<class Type>
+Type maxMagSqr(const UList<Type>& f)
+{
+    if (f.size())
+    {
+        Type Max(f[0]);
+        TFOR_ALL_S_OP_FUNC_F_S
+        (
+            Type,
+            Max,
+            =,
+            maxMagSqrOp<Type>(),
+            Type,
+            f,
+            Type,
+            Max
+        )
+        return Max;
+    }
+    else
+    {
+        return pTraits<Type>::min;
+    }
+}
+
+TMP_UNARY_FUNCTION(Type, maxMagSqr)
+
+template<class Type>
+Type minMagSqr(const UList<Type>& f)
+{
+    if (f.size())
+    {
+        Type Min(f[0]);
+        TFOR_ALL_S_OP_FUNC_F_S
+        (
+            Type,
+            Min,
+            =,
+            minMagSqrOp<Type>(),
+            Type,
+            f,
+            Type,
+            Min
+        )
+        return Min;
+    }
+    else
+    {
+        return pTraits<Type>::max;
+    }
+}
+
+TMP_UNARY_FUNCTION(Type, minMagSqr)
 
 template<class Type>
 scalar sumProd(const UList<Type>& f1, const UList<Type>& f2)
@@ -488,6 +541,8 @@ TMP_UNARY_FUNCTION(ReturnType, gFunc)
 G_UNARY_FUNCTION(Type, gMax, max, max)
 G_UNARY_FUNCTION(Type, gMin, min, min)
 G_UNARY_FUNCTION(Type, gSum, sum, sum)
+G_UNARY_FUNCTION(Type, gMaxMagSqr, maxMagSqr, maxMagSqr)
+G_UNARY_FUNCTION(Type, gMinMagSqr, minMagSqr, minMagSqr)
 G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
 G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
 G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
index f1fa2fdd968281bd575b0c40152099d110078705..4c51262f9e36c9f76e1489c1ab31cb6280ca0348 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
+++ b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
@@ -171,6 +171,16 @@ Type sum(const UList<Type>& f);
 
 TMP_UNARY_FUNCTION(Type, sum)
 
+template<class Type>
+Type maxMagSqr(const UList<Type>& f);
+
+TMP_UNARY_FUNCTION(Type, maxMagSqr)
+
+template<class Type>
+Type minMagSqr(const UList<Type>& f);
+
+TMP_UNARY_FUNCTION(Type, minMagSqr)
+
 
 template<class Type>
 scalar sumProd(const UList<Type>& f1, const UList<Type>& f2);
@@ -208,6 +218,8 @@ TMP_UNARY_FUNCTION(ReturnType, gFunc)
 G_UNARY_FUNCTION(Type, gMax, max, max)
 G_UNARY_FUNCTION(Type, gMin, min, min)
 G_UNARY_FUNCTION(Type, gSum, sum, sum)
+G_UNARY_FUNCTION(Type, gMaxMagSqr, maxMagSqr, maxMagSqr)
+G_UNARY_FUNCTION(Type, gMinMagSqr, minMagSqr, minMagSqr)
 G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
 G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
 G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index a83b2bf1909f9e962ac843bd95e981d290148a92..2739e828b349f34986ea6028e57f2967de45f885 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -1963,7 +1963,7 @@ Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
     pointField sharedPoints(mesh_.points(), sharedPointLabels());
 
     // Append from all processors
-    combineReduce(sharedPoints, plusEqOp<pointField>());
+    combineReduce(sharedPoints, ListPlusEqOp<pointField>());
 
     // Merge tolerance
     scalar tolDim = matchTol_ * mesh_.bounds().mag();
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
index ab1858b1484f436d7434f0e5c84044cb34a15b9e..35543fd0937d399aa4df018d4f0e75b88c0c8482 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
@@ -109,29 +109,6 @@ class globalMeshData
     public processorTopology
 {
 
-    // Private class
-
-        // To combineReduce a pointField. Just appends all lists.
-        template<class T>
-        class plusEqOp
-        {
-
-        public:
-
-            void operator()(T& x, const T& y) const
-            {
-                label n = x.size();
-
-                x.setSize(x.size() + y.size());
-
-                forAll(y, i)
-                {
-                    x[n++] = y[i];
-                }
-            }
-        };
-
-
     // Private data
 
         //- Reference to mesh
@@ -337,6 +314,29 @@ class globalMeshData
 
 public:
 
+    // Public class
+
+        // To combineReduce a List. Just appends all lists.
+        template<class T>
+        class ListPlusEqOp
+        {
+
+        public:
+
+            void operator()(T& x, const T& y) const
+            {
+                label n = x.size();
+
+                x.setSize(x.size() + y.size());
+
+                forAll(y, i)
+                {
+                    x[n++] = y[i];
+                }
+            }
+        };
+
+
     //- Runtime type information
     ClassName("globalMeshData");
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 685e043fdb02e08740120223cdfa50bcce64432e..a0a325410d2b70fe7de814fac0f9099f4477a537 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -1446,9 +1446,9 @@ void Foam::autoLayerDriver::syncPatchDisplacement
         }
     }
 
-    Info<< "Prevented extrusion on "
-        << returnReduce(nChangedTotal, sumOp<label>())
-        << " coupled patch points during syncPatchDisplacement." << endl;
+    //Info<< "Prevented extrusion on "
+    //    << returnReduce(nChangedTotal, sumOp<label>())
+    //    << " coupled patch points during syncPatchDisplacement." << endl;
 }
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
index 9ce03ec25c80b646eefc4971f6f967b55ea5b25f..3e824b1e0287624034d6af5e6e993783a6164e2b 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
@@ -524,7 +524,7 @@ void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
 
 
 
-    label nOldPointCounter = nPointCounter;
+    //label nOldPointCounter = nPointCounter;
 
     // Detect situation where two featureedge-neighbouring faces are partly or
     // not extruded and the edge itself is extruded. In this case unmark the
@@ -633,9 +633,9 @@ void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
         }
     }
 
-    Info<< "Added "
-        << returnReduce(nPointCounter-nOldPointCounter, sumOp<label>())
-        << " point not to extrude." << endl;
+    //Info<< "Added "
+    //    << returnReduce(nPointCounter-nOldPointCounter, sumOp<label>())
+    //    << " point not to extrude." << endl;
 }
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 07ff697df6a010fed0e6e41bcd3eeccf3f81f971..adb7c3859032a928765510d046141f802df54794 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -910,6 +910,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
     Info<< "Detecting near surfaces ..." << endl;
 
     const pointField& localPoints = pp.localPoints();
+    const labelList& meshPoints = pp.meshPoints();
     const refinementSurfaces& surfaces = meshRefiner_.surfaces();
     const fvMesh& mesh = meshRefiner_.mesh();
 
@@ -1220,6 +1221,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
     }
 
 
+    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
     label nOverride = 0;
 
     // 1. All points to non-interface surfaces
@@ -1332,7 +1334,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
                 }
             }
 
-            if (override)
+            if (override && isMasterPoint[meshPoints[pointI]])
             {
                 nOverride++;
             }
@@ -1399,8 +1401,6 @@ void Foam::autoSnapDriver::detectNearSurfaces
             );
 
 
-            label nOverride = 0;
-
             forAll(hit1, i)
             {
                 label pointI = zonePointIndices[i];
@@ -1472,7 +1472,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
                     }
                 }
 
-                if (override)
+                if (override && isMasterPoint[meshPoints[pointI]])
                 {
                     nOverride++;
                 }
@@ -1690,7 +1690,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             scalarField magDisp(mag(patchDisp));
 
             Info<< "Wanted displacement : average:"
-                << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
+                <<  meshRefinement::gAverage
+                    (
+                        mesh,
+                        syncTools::getMasterPoints(mesh),
+                        pp.meshPoints(),
+                        magDisp
+                    )
                 << " min:" << gMin(magDisp)
                 << " max:" << gMax(magDisp) << endl;
         }
@@ -2548,25 +2554,30 @@ void Foam::autoSnapDriver::doSnap
                 adaptPatchIDs
             )
         );
-        indirectPrimitivePatch& pp = ppPtr();
+
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp));
+        const scalarField snapDist(calcSnapDistance(mesh, snapParams, ppPtr()));
 
 
         // Construct iterative mesh mover.
         Info<< "Constructing mesh displacer ..." << endl;
         Info<< "Using mesh parameters " << motionDict << nl << endl;
 
-        const pointMesh& pMesh = pointMesh::New(mesh);
-
-        motionSmoother meshMover
+        autoPtr<motionSmoother> meshMoverPtr
         (
-            mesh,
-            pp,
-            adaptPatchIDs,
-            meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
-            motionDict
+            new motionSmoother
+            (
+                mesh,
+                ppPtr(),
+                adaptPatchIDs,
+                meshRefinement::makeDisplacementField
+                (
+                    pointMesh::New(mesh),
+                    adaptPatchIDs
+                ),
+                motionDict
+            )
         );
 
 
@@ -2595,16 +2606,95 @@ void Foam::autoSnapDriver::doSnap
             snapParams,
             nInitErrors,
             baffles,
-            meshMover
+            meshMoverPtr()
         );
 
 
+
+        //- Only if in feature attraction mode:
+        // Nearest feature
+        vectorField patchAttraction;
+        // Constraints at feature
+        List<pointConstraint> patchConstraints;
+
+
         for (label iter = 0; iter < nFeatIter; iter++)
         {
             Info<< nl
                 << "Morph iteration " << iter << nl
                 << "-----------------" << endl;
 
+
+            //if (iter > 0 && iter == nFeatIter/2)
+            //{
+            //    Info<< "Splitting diagonal attractions" << endl;
+            //    const labelList& bFaces = ppPtr().addressing();
+            //
+            //    DynamicList<label> splitFaces(bFaces.size());
+            //    DynamicList<labelPair> splits(bFaces.size());
+            //
+            //    forAll(bFaces, faceI)
+            //    {
+            //        const labelPair split
+            //        (
+            //            findDiagonalAttraction
+            //            (
+            //                ppPtr(),
+            //                patchAttraction,
+            //                patchConstraints,
+            //                faceI
+            //            )
+            //        );
+            //
+            //        if (split != labelPair(-1, -1))
+            //        {
+            //            splitFaces.append(bFaces[faceI]);
+            //            splits.append(split);
+            //        }
+            //    }
+            //
+            //    autoPtr<mapPolyMesh> mapPtr = meshRefiner_.splitFaces
+            //    (
+            //        splitFaces,
+            //        splits
+            //    );
+            //
+            //    const labelList& faceMap = mapPtr().faceMap();
+            //    meshRefinement::updateList(faceMap, -1, duplicateFace);
+            //    const labelList& reverseFaceMap = mapPtr().reverseFaceMap();
+            //    forAll(baffles, i)
+            //    {
+            //        labelPair& baffle = baffles[i];
+            //        baffle.first() = reverseFaceMap[baffle.first()];
+            //        baffle.second() = reverseFaceMap[baffle.second()];
+            //    }
+            //
+            //    meshMoverPtr.clear();
+            //    ppPtr.clear();
+            //
+            //    ppPtr = meshRefinement::makePatch(mesh, adaptPatchIDs);
+            //    meshMoverPtr.reset
+            //    (
+            //        new motionSmoother
+            //        (
+            //            mesh,
+            //            ppPtr(),
+            //            adaptPatchIDs,
+            //            meshRefinement::makeDisplacementField
+            //            (
+            //                pointMesh::New(mesh),
+            //                adaptPatchIDs
+            //            ),
+            //            motionDict
+            //        )
+            //    );
+            //}
+
+
+            indirectPrimitivePatch& pp = ppPtr();
+            motionSmoother& meshMover = meshMoverPtr();
+
+
             // Calculate displacement at every patch point. Insert into
             // meshMover.
             // Calculate displacement at every patch point
@@ -2652,7 +2742,9 @@ void Foam::autoSnapDriver::doSnap
                     scalar(iter+1)/nFeatIter,
                     snapDist,
                     disp,
-                    meshMover
+                    meshMover,
+                    patchAttraction,
+                    patchConstraints
                 );
             }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 77327c0689b0728b34f4201c0a09af4a9d874650..b8282bedfcd8e5f0399be522e5eccbfe2b1bf827 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -125,7 +125,9 @@ class autoSnapDriver
 
                 void smoothAndConstrain
                 (
+                    const PackedBoolList& isMasterEdge,
                     const indirectPrimitivePatch& pp,
+                    const labelList& meshEdges,
                     const List<pointConstraint>& constraints,
                     vectorField& disp
                 ) const;
@@ -196,6 +198,16 @@ class autoSnapDriver
                     List<pointConstraint>& patchConstraints
                 ) const;
 
+                //- Detect any diagonal attraction. Returns indices in face
+                //  or (-1, -1) if none
+                labelPair findDiagonalAttraction
+                (
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& patchAttraction,
+                    const List<pointConstraint>& patchConstraints,
+                    const label faceI
+                ) const;
+
                 //- Return hit if on multiple points
                 pointIndexHit findMultiPatchPoint
                 (
@@ -371,7 +383,9 @@ class autoSnapDriver
                     const scalar featureAttract,
                     const scalarField& snapDist,
                     const vectorField& nearestDisp,
-                    motionSmoother& meshMover
+                    motionSmoother& meshMover,
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
                 ) const;
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 887674ea7aebe982742e249b2681c1791d8b9b50..32133ad5f8437c3e754777792bf6385a1260c951 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -38,50 +38,18 @@ License
 #include "treeDataPoint.H"
 #include "indexedOctree.H"
 #include "snapParameters.H"
+#include "PatchTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    class listTransform
-    {
-    public:
-
-        void operator()
-        (
-            const vectorTensorTransform& vt,
-            const bool forward,
-            List<List<point> >& fld
-        ) const
-        {
-            const tensor T
-            (
-                forward
-              ? vt.R()
-              : vt.R().T()
-            );
-
-            forAll(fld, i)
-            {
-                List<point>& elems = fld[i];
-                forAll(elems, elemI)
-                {
-                    elems[elemI] = transform(T, elems[elemI]);
-                }
-            }
-        }
-    };
-
     template<class T>
     class listPlusEqOp
     {
     public:
 
-        void operator()
-        (
-            List<T>& x,
-            const List<T>& y
-        ) const
+        void operator()(List<T>& x, const List<T>& y) const
         {
             label sz = x.size();
             x.setSize(sz+y.size());
@@ -159,7 +127,9 @@ bool Foam::autoSnapDriver::isFeaturePoint
 
 void Foam::autoSnapDriver::smoothAndConstrain
 (
+    const PackedBoolList& isMasterEdge,
     const indirectPrimitivePatch& pp,
+    const labelList& meshEdges,
     const List<pointConstraint>& constraints,
     vectorField& disp
 ) const
@@ -197,11 +167,16 @@ void Foam::autoSnapDriver::smoothAndConstrain
             {
                 forAll(pEdges, i)
                 {
-                    label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
-                    if (constraints[nbrPointI].first() >= nConstraints)
+                    label edgeI = pEdges[i];
+
+                    if (isMasterEdge[meshEdges[edgeI]])
                     {
-                        dispSum[pointI] += disp[nbrPointI];
-                        dispCount[pointI]++;
+                        label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
+                        if (constraints[nbrPointI].first() >= nConstraints)
+                        {
+                            dispSum[pointI] += disp[nbrPointI];
+                            dispCount[pointI]++;
+                        }
                     }
                 }
             }
@@ -564,6 +539,9 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
 {
     const fvMesh& mesh = meshRefiner_.mesh();
 
+    const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
+
+
     // For now just get all surrounding face data. Expensive - should just
     // store and sync data on coupled points only
     // (see e.g PatchToolsNormals.C)
@@ -583,7 +561,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         forAll(pFaces, i)
         {
             label faceI = pFaces[i];
-            if (faceSurfaceGlobalRegion[faceI] != -1)
+            if (isMasterFace[faceI] && faceSurfaceGlobalRegion[faceI] != -1)
             {
                 nFaces++;
             }
@@ -605,7 +583,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
             label faceI = pFaces[i];
             label globalRegionI = faceSurfaceGlobalRegion[faceI];
 
-            if (globalRegionI != -1)
+            if (isMasterFace[faceI] && globalRegionI != -1)
             {
                 pNormals[nFaces] = faceSurfaceNormal[faceI];
                 pDisp[nFaces] = faceDisp[faceI];
@@ -694,7 +672,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pointFaceSurfNormals,
         listPlusEqOp<point>(),
         List<point>(),
-        listTransform()
+        mapDistribute::transform()
     );
     syncTools::syncPointList
     (
@@ -703,7 +681,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pointFaceDisp,
         listPlusEqOp<point>(),
         List<point>(),
-        listTransform()
+        mapDistribute::transform()
     );
     syncTools::syncPointList
     (
@@ -712,7 +690,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pointFaceCentres,
         listPlusEqOp<point>(),
         List<point>(),
-        listTransform()
+        mapDistribute::transformPosition()
     );
     syncTools::syncPointList
     (
@@ -722,6 +700,25 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         listPlusEqOp<label>(),
         List<label>()
     );
+
+
+    // Sort the data according to the face centres. This is only so we get
+    // consistent behaviour serial and parallel.
+    labelList visitOrder;
+    forAll(pointFaceDisp, pointI)
+    {
+        List<point>& pNormals = pointFaceSurfNormals[pointI];
+        List<point>& pDisp = pointFaceDisp[pointI];
+        List<point>& pFc = pointFaceCentres[pointI];
+        labelList& pFid = pointFacePatchID[pointI];
+
+        sortedOrder(mag(pFc)(), visitOrder);
+
+        pNormals = List<point>(pNormals, visitOrder);
+        pDisp = List<point>(pDisp, visitOrder);
+        pFc = List<point>(pFc, visitOrder);
+        pFid = UIndirectList<label>(pFid, visitOrder);
+    }
 }
 
 
@@ -1360,6 +1357,70 @@ void Foam::autoSnapDriver::stringFeatureEdges
 }
 
 
+// If only two attractions and across face return the face indices
+Foam::labelPair Foam::autoSnapDriver::findDiagonalAttraction
+(
+    const indirectPrimitivePatch& pp,
+    const vectorField& patchAttraction,
+    const List<pointConstraint>& patchConstraints,
+    const label faceI
+) const
+{
+    const face& f = pp.localFaces()[faceI];
+    // For now just detect any attraction. Improve this to look at
+    // actual attraction position and orientation
+
+    labelPair attractIndices(-1, -1);
+
+    if (f.size() >= 4)
+    {
+        forAll(f, fp)
+        {
+            label pointI = f[fp];
+            if (patchConstraints[pointI].first() >= 2)
+            {
+                // Attract to feature edge or point
+                if (attractIndices[0] == -1)
+                {
+                    // First attraction. Store
+                    attractIndices[0] = fp;
+                }
+                else if (attractIndices[1] == -1)
+                {
+                    // Second attraction. Check if not consecutive to first
+                    // attraction
+                    label fp0 = attractIndices[0];
+                    if (f.fcIndex(fp0) == fp || f.fcIndex(fp) == fp0)
+                    {
+                        // Consecutive. Skip.
+                        attractIndices = labelPair(-1, -1);
+                        break;
+                    }
+                    else
+                    {
+                        attractIndices[1] = fp;
+                    }
+                }
+                else
+                {
+                    // More than two attractions. Skip.
+                    attractIndices = labelPair(-1, -1);
+                    break;
+                }
+            }
+        }
+
+
+        if (attractIndices[1] == -1)
+        {
+            // Found only one attraction. Skip.
+            attractIndices = labelPair(-1, -1);
+        }
+    }
+    return attractIndices;
+}
+
+
 Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
 (
     const indirectPrimitivePatch& pp,
@@ -1929,7 +1990,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                 edgeFaceNormals,
                 listPlusEqOp<point>(),
                 List<point>(),
-                listTransform()
+                mapDistribute::transform()
             );
         }
 
@@ -2778,7 +2839,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     const scalar featureAttract,
     const scalarField& snapDist,
     const vectorField& nearestDisp,
-    motionSmoother& meshMover
+    motionSmoother& meshMover,
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
 ) const
 {
     const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
@@ -2851,7 +2914,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // - faceSurfaceNormal
     // - faceDisp
-    // - faceCentres&faceNormal
+    // - faceCentres
     List<List<point> > pointFaceSurfNormals(pp.nPoints());
     List<List<point> > pointFaceDisp(pp.nPoints());
     List<List<point> > pointFaceCentres(pp.nPoints());
@@ -2884,10 +2947,11 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // here.
 
     // Nearest feature
-    vectorField patchAttraction(localPoints.size(), vector::zero);
+    patchAttraction.setSize(localPoints.size());
+    patchAttraction = vector::zero;
     // Constraints at feature
-    List<pointConstraint> patchConstraints(localPoints.size());
-
+    patchConstraints.setSize(localPoints.size());
+    patchConstraints = pointConstraint();
 
     if (implicitFeatureAttraction)
     {
@@ -2951,15 +3015,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         patchConstraints
     );
 
+    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+    {
+        vector avgPatchDisp = meshRefinement::gAverage
+        (
+            mesh,
+            isMasterPoint,
+            pp.meshPoints(),
+            patchDisp
+        );
+        vector avgPatchAttr = meshRefinement::gAverage
+        (
+            mesh,
+            isMasterPoint,
+            pp.meshPoints(),
+            patchAttraction
+        );
 
-    Info<< "Attraction:" << endl
-        << "     linear   : max:" << gMax(patchDisp)
-        << " avg:" << gAverage(patchDisp)
-        << endl
-        << "     feature  : max:" << gMax(patchAttraction)
-        << " avg:" << gAverage(patchAttraction)
-        << endl;
-
+        Info<< "Attraction:" << endl
+            << "     linear   : max:" << gMaxMagSqr(patchDisp)
+            << " avg:" << avgPatchDisp << endl
+            << "     feature  : max:" << gMaxMagSqr(patchAttraction)
+            << " avg:" << avgPatchAttr << endl;
+    }
 
     // So now we have:
     // - patchDisp          : point movement to go to nearest point on surface
@@ -2986,37 +3064,45 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     // Count
     {
+        const labelList& meshPoints = pp.meshPoints();
+
+        label nMasterPoints = 0;
         label nPlanar = 0;
         label nEdge = 0;
         label nPoint = 0;
 
         forAll(patchConstraints, pointI)
         {
-            if (patchConstraints[pointI].first() == 1)
+            if (isMasterPoint[meshPoints[pointI]])
             {
-                nPlanar++;
-            }
-            else if (patchConstraints[pointI].first() == 2)
-            {
-                nEdge++;
-            }
-            else if (patchConstraints[pointI].first() == 3)
-            {
-                nPoint++;
+                nMasterPoints++;
+
+                if (patchConstraints[pointI].first() == 1)
+                {
+                    nPlanar++;
+                }
+                else if (patchConstraints[pointI].first() == 2)
+                {
+                    nEdge++;
+                }
+                else if (patchConstraints[pointI].first() == 3)
+                {
+                    nPoint++;
+                }
             }
         }
 
-        label nTotPoints = returnReduce(pp.nPoints(), sumOp<label>());
+        reduce(nMasterPoints, sumOp<label>());
         reduce(nPlanar, sumOp<label>());
         reduce(nEdge, sumOp<label>());
         reduce(nPoint, sumOp<label>());
-        Info<< "Feature analysis : total points:"
-            << nTotPoints
+        Info<< "Feature analysis : total master points:"
+            << nMasterPoints
             << " attraction to :" << nl
             << "    feature point   : " << nPoint << nl
             << "    feature edge    : " << nEdge << nl
             << "    nearest surface : " << nPlanar << nl
-            << "    rest            : " << nTotPoints-nPoint-nEdge-nPlanar
+            << "    rest            : " << nMasterPoints-nPoint-nEdge-nPlanar
             << nl
             << endl;
     }
@@ -3035,11 +3121,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     if (featureAttract < 1-0.001)
     {
+        const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
+
+        const vectorField pointNormals
+        (
+            PatchTools::pointNormals
+            (
+                mesh,
+                pp
+            )
+        );
+        const labelList meshEdges
+        (
+            pp.meshEdges(mesh.edges(), mesh.pointEdges())
+        );
+
+
         // 1. Smoothed all displacement
         vectorField smoothedPatchDisp = patchDisp;
         smoothAndConstrain
         (
+            isMasterEdge,
             pp,
+            meshEdges,
             patchConstraints,
             smoothedPatchDisp
         );
@@ -3047,16 +3151,18 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
         // 2. Smoothed tangential component
         vectorField tangPatchDisp = patchDisp;
-        tangPatchDisp -= (pp.pointNormals() & patchDisp) * pp.pointNormals();
+        tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
         smoothAndConstrain
         (
+            isMasterEdge,
             pp,
+            meshEdges,
             patchConstraints,
             tangPatchDisp
         );
 
         // Re-add normal component
-        tangPatchDisp += (pp.pointNormals() & patchDisp) * pp.pointNormals();
+        tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
 
         if (debug&meshRefinement::OBJINTERSECTIONS)
         {
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 8a5f327b4108229c61ca28ea5a177ecc48c7f356..873167bc4f0227cfb7302e7f79ef7382c6e3e8cb 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -259,23 +259,14 @@ private:
                 label& nRefine
             ) const;
 
-            //- Mark cell if local curvature > curvature or
-            //  markDifferingRegions = true and intersections with different
-            //  regions.
-            bool checkCurvature
+            //- Helper: count number of normals1 that are in normals2
+            label countMatches
             (
-                const scalar curvature,
-                const label nAllowRefine,
-                const label surfaceLevel,
-                const vector& surfaceNormal,
-                const label cellI,
-                label& cellMaxLevel,
-                vector& cellMaxNormal,
-                labelList& refineCell,
-                label& nRefine
+                const List<point>& normals1,
+                const List<point>& normals2,
+                const scalar tol = 1e-6
             ) const;
 
-
             //- Mark cells for surface curvature based refinement. Marks if
             //  local curvature > curvature or if on different regions
             //  (markDifferingRegions)
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 7444493e2982bfe088d20f804db3facc52b5c3b9..ce84ceec38eeb6c633986d83de12bb7a2a1fdb4d 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -38,7 +38,74 @@ License
 #include "featureEdgeMesh.H"
 #include "Cloud.H"
 //#include "globalIndex.H"
-//#include "OBJstream.H"
+#include "OBJstream.H"
+#include "cellSet.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    //- To compare normals
+    static bool less(const vector& x, const vector& y)
+    {
+        for (direction i = 0; i < vector::nComponents; i++)
+        {
+            if (x[i] < y[i])
+            {
+                return true;
+            }
+            else if (x[i] > y[i])
+            {
+                return false;
+            }
+        }
+        // All components the same
+        return false;
+    }
+
+
+    //- To compare normals
+    class normalLess
+    {
+        const vectorList& values_;
+
+    public:
+
+        normalLess(const vectorList& values)
+        :
+            values_(values)
+        {}
+
+        bool operator()(const label a, const label b) const
+        {
+            return less(values_[a], values_[b]);
+        }
+    };
+
+
+    //- template specialization for pTraits<labelList> so we can have fields
+    template<>
+    class pTraits<labelList>
+    {
+
+    public:
+
+        //- Component type
+        typedef labelList cmptType;
+    };
+
+    //- template specialization for pTraits<labelList> so we can have fields
+    template<>
+    class pTraits<vectorList>
+    {
+
+    public:
+
+        //- Component type
+        typedef vectorList cmptType;
+    };
+}
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -954,64 +1021,33 @@ Foam::label Foam::meshRefinement::markSurfaceRefinement
 }
 
 
-// Checks if multiple intersections of a cell (by a surface with a higher
-// max than the cell level) and if so if the normals at these intersections
-// make a large angle.
-// Returns false if the nRefine limit has been reached, true otherwise.
-bool Foam::meshRefinement::checkCurvature
+// Count number of matches of first argument in second argument
+Foam::label Foam::meshRefinement::countMatches
 (
-    const scalar curvature,
-    const label nAllowRefine,
-
-    const label surfaceLevel,   // current intersection max level
-    const vector& surfaceNormal,// current intersection normal
-
-    const label cellI,
-
-    label& cellMaxLevel,        // cached max surface level for this cell
-    vector& cellMaxNormal,      // cached surface normal for this cell
-
-    labelList& refineCell,
-    label& nRefine
+    const List<point>& normals1,
+    const List<point>& normals2,
+    const scalar tol
 ) const
 {
-    const labelList& cellLevel = meshCutter_.cellLevel();
+    label nMatches = 0;
 
-    // Test if surface applicable
-    if (surfaceLevel > cellLevel[cellI])
+    forAll(normals1, i)
     {
-        if (cellMaxLevel == -1)
-        {
-            // First visit of cell. Store
-            cellMaxLevel = surfaceLevel;
-            cellMaxNormal = surfaceNormal;
-        }
-        else
+        const vector& n1 = normals1[i];
+
+        forAll(normals2, j)
         {
-            // Second or more visit. Check curvature.
-            if ((cellMaxNormal & surfaceNormal) < curvature)
-            {
-                return markForRefine
-                (
-                    surfaceLevel,   // mark with any non-neg number.
-                    nAllowRefine,
-                    refineCell[cellI],
-                    nRefine
-                );
-            }
+            const vector& n2 = normals2[j];
 
-            // Set normal to that of highest surface. Not really necessary
-            // over here but we reuse cellMax info when doing coupled faces.
-            if (surfaceLevel > cellMaxLevel)
+            if (magSqr(n1-n2) < tol)
             {
-                cellMaxLevel = surfaceLevel;
-                cellMaxNormal = surfaceNormal;
+                nMatches++;
+                break;
             }
         }
     }
 
-    // Did not reach refinement limit.
-    return true;
+    return nMatches;
 }
 
 
@@ -1039,6 +1075,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
     // on a different surface gets refined (if its current level etc.)
 
 
+    const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
+
+
     // Collect candidate faces (i.e. intersecting any surface and
     // owner/neighbour not yet refined.
     labelList testFaces(getRefineCandidateFaces(refineCell));
@@ -1068,6 +1107,12 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 
             start[i] = cellCentres[own];
             end[i] = neiCc[bFaceI];
+
+            if (!isMasterFace[faceI])
+            {
+                Swap(start[i], end[i]);
+            }
+
             minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
         }
     }
@@ -1081,10 +1126,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 
 
     // Test for all intersections (with surfaces of higher max level than
-    // minLevel) and cache per cell the max surface level and the local normal
-    // on that surface.
-    labelList cellMaxLevel(mesh_.nCells(), -1);
-    vectorField cellMaxNormal(mesh_.nCells(), vector::zero);
+    // minLevel) and cache per cell the interesting inter
+    labelListList cellSurfLevels(mesh_.nCells());
+    List<vectorList> cellSurfNormals(mesh_.nCells());
 
     {
         // Per segment the normals of the surfaces
@@ -1104,12 +1148,29 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
             surfaceNormal,
             surfaceLevel
         );
+
+
+        // Sort the data according to surface location. This will guarantee
+        // that on coupled faces both sides visit the intersections in
+        // the same order so will decide the same
+        labelList visitOrder;
+        forAll(surfaceNormal, pointI)
+        {
+            vectorList& pNormals = surfaceNormal[pointI];
+            labelList& pLevel = surfaceLevel[pointI];
+
+            sortedOrder(pNormals, visitOrder, normalLess(pNormals));
+
+            pNormals = List<point>(pNormals, visitOrder);
+            pLevel = UIndirectList<label>(pLevel, visitOrder);
+        }
+
         // Clear out unnecessary data
         start.clear();
         end.clear();
         minLevel.clear();
 
-        // Extract per cell information on the surface with the highest max
+        // Convert face-wise data to cell.
         forAll(testFaces, i)
         {
             label faceI = testFaces[i];
@@ -1118,164 +1179,280 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
             const vectorList& fNormals = surfaceNormal[i];
             const labelList& fLevels = surfaceLevel[i];
 
-            forAll(fLevels, hitI)
+            forAll(fNormals, hitI)
             {
-                checkCurvature
-                (
-                    curvature,
-                    nAllowRefine,
-
-                    fLevels[hitI],
-                    fNormals[hitI],
-
-                    own,
-                    cellMaxLevel[own],
-                    cellMaxNormal[own],
+                if (fLevels[hitI] > cellLevel[own])
+                {
+                    cellSurfLevels[own].append(fLevels[hitI]);
+                    cellSurfNormals[own].append(fNormals[hitI]);
+                }
 
-                    refineCell,
-                    nRefine
-                );
+                if (mesh_.isInternalFace(faceI))
+                {
+                    label nei = mesh_.faceNeighbour()[faceI];
+                    if (fLevels[hitI] > cellLevel[nei])
+                    {
+                        cellSurfLevels[nei].append(fLevels[hitI]);
+                        cellSurfNormals[nei].append(fNormals[hitI]);
+                    }
+                }
             }
+        }
+    }
 
-            if (mesh_.isInternalFace(faceI))
-            {
-                label nei = mesh_.faceNeighbour()[faceI];
-
-                forAll(fLevels, hitI)
-                {
-                    checkCurvature
-                    (
-                        curvature,
-                        nAllowRefine,
 
-                        fLevels[hitI],
-                        fNormals[hitI],
 
-                        nei,
-                        cellMaxLevel[nei],
-                        cellMaxNormal[nei],
-
-                        refineCell,
-                        nRefine
-                    );
-                }
+    // Bit of statistics
+    {
+        label nSet = 0;
+        label nNormals = 9;
+        forAll(cellSurfNormals, cellI)
+        {
+            const vectorList& normals = cellSurfNormals[cellI];
+            if (normals.size())
+            {
+                nSet++;
+                nNormals += normals.size();
             }
         }
+        reduce(nSet, sumOp<label>());
+        reduce(nNormals, sumOp<label>());
+        Info<< "markSurfaceCurvatureRefinement :"
+            << " cells:" << mesh_.globalData().nTotalCells()
+            << " of which with normals:" << nSet
+            << " ; total normals stored:" << nNormals
+            << endl;
     }
 
-    // 2. Find out a measure of surface curvature and region edges.
-    // Send over surface region and surface normal to neighbour cell.
 
-    labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
 
-    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
-    {
-        label bFaceI = faceI-mesh_.nInternalFaces();
-        label own = mesh_.faceOwner()[faceI];
+    bool reachedLimit = false;
 
-        neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
-        neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
-    }
-    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
-    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
 
-    // Loop over all faces. Could only be checkFaces.. except if they're coupled
+    // 1. Check normals per cell
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Internal faces
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    for
+    (
+        label cellI = 0;
+        !reachedLimit && cellI < cellSurfNormals.size();
+        cellI++
+    )
     {
-        label own = mesh_.faceOwner()[faceI];
-        label nei = mesh_.faceNeighbour()[faceI];
+        const vectorList& normals = cellSurfNormals[cellI];
+        const labelList& levels = cellSurfLevels[cellI];
 
-        if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
+        // n^2 comparison of all normals in a cell
+        for (label i = 0; !reachedLimit && i < normals.size(); i++)
         {
-            // Have valid data on both sides. Check curvature.
-            if ((cellMaxNormal[own] & cellMaxNormal[nei]) < curvature)
+            for (label j = i+1; !reachedLimit && j < normals.size(); j++)
             {
-                // See which side to refine
-                if (cellLevel[own] < cellMaxLevel[own])
+                if ((normals[i] & normals[j]) < curvature)
                 {
-                    if
-                    (
-                        !markForRefine
+                    label maxLevel = max(levels[i], levels[j]);
+
+                    if (cellLevel[cellI] < maxLevel)
+                    {
+                        if
                         (
-                            cellMaxLevel[own],
-                            nAllowRefine,
-                            refineCell[own],
-                            nRefine
+                            !markForRefine
+                            (
+                                maxLevel,
+                                nAllowRefine,
+                                refineCell[cellI],
+                                nRefine
+                            )
                         )
-                    )
-                    {
-                        if (debug)
                         {
-                            Pout<< "Stopped refining since reaching my cell"
-                                << " limit of " << mesh_.nCells()+7*nRefine
-                                << endl;
+                            if (debug)
+                            {
+                                Pout<< "Stopped refining since reaching my cell"
+                                    << " limit of " << mesh_.nCells()+7*nRefine
+                                    << endl;
+                            }
+                            reachedLimit = true;
+                            break;
                         }
-                        break;
                     }
                 }
+            }
+        }
+    }
 
-                if (cellLevel[nei] < cellMaxLevel[nei])
+
+
+    // 2. Find out a measure of surface curvature
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Look at normals between neighbouring surfaces
+    // Loop over all faces. Could only be checkFaces, except if they're coupled
+
+    // Internal faces
+    for
+    (
+        label faceI = 0;
+        !reachedLimit && faceI < mesh_.nInternalFaces();
+        faceI++
+    )
+    {
+        label own = mesh_.faceOwner()[faceI];
+        label nei = mesh_.faceNeighbour()[faceI];
+
+        const vectorList& ownNormals = cellSurfNormals[own];
+        const labelList& ownLevels = cellSurfLevels[own];
+        const vectorList& neiNormals = cellSurfNormals[nei];
+        const labelList& neiLevels = cellSurfLevels[nei];
+
+
+        // Special case: owner normals set is a subset of the neighbour
+        // normals. Do not do curvature refinement since other cell's normals
+        // do not add any information. Avoids weird corner extensions of
+        // refinement regions.
+        bool ownIsSubset =
+            countMatches(ownNormals, neiNormals)
+         == ownNormals.size();
+
+        bool neiIsSubset =
+            countMatches(neiNormals, ownNormals)
+         == neiNormals.size();
+
+
+        if (!ownIsSubset && !neiIsSubset)
+        {
+            // n^2 comparison of between ownNormals and neiNormals
+            for (label i = 0; !reachedLimit &&  i < ownNormals.size(); i++)
+            {
+                for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
                 {
-                    if
-                    (
-                        !markForRefine
-                        (
-                            cellMaxLevel[nei],
-                            nAllowRefine,
-                            refineCell[nei],
-                            nRefine
-                        )
-                    )
+                    // Have valid data on both sides. Check curvature.
+                    if ((ownNormals[i] & neiNormals[j]) < curvature)
                     {
-                        if (debug)
+                        // See which side to refine.
+                        if (cellLevel[own] < ownLevels[i])
                         {
-                            Pout<< "Stopped refining since reaching my cell"
-                                << " limit of " << mesh_.nCells()+7*nRefine
-                                << endl;
+                            if
+                            (
+                                !markForRefine
+                                (
+                                    ownLevels[i],
+                                    nAllowRefine,
+                                    refineCell[own],
+                                    nRefine
+                                )
+                            )
+                            {
+                                if (debug)
+                                {
+                                    Pout<< "Stopped refining since reaching"
+                                        << " my cell limit of "
+                                        << mesh_.nCells()+7*nRefine << endl;
+                                }
+                                reachedLimit = true;
+                                break;
+                            }
+                        }
+                        if (cellLevel[nei] < neiLevels[j])
+                        {
+                            if
+                            (
+                                !markForRefine
+                                (
+                                    neiLevels[j],
+                                    nAllowRefine,
+                                    refineCell[nei],
+                                    nRefine
+                                )
+                            )
+                            {
+                                if (debug)
+                                {
+                                    Pout<< "Stopped refining since reaching"
+                                        << " my cell limit of "
+                                        << mesh_.nCells()+7*nRefine << endl;
+                                }
+                                reachedLimit = true;
+                                break;
+                            }
                         }
-                        break;
                     }
                 }
             }
         }
     }
+
+
+    // Send over surface normal to neighbour cell.
+    List<vectorList> neiSurfaceNormals;
+    syncTools::swapBoundaryCellList(mesh_, cellSurfNormals, neiSurfaceNormals);
+
     // Boundary faces
-    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+    for
+    (
+        label faceI = mesh_.nInternalFaces();
+        !reachedLimit && faceI < mesh_.nFaces();
+        faceI++
+    )
     {
         label own = mesh_.faceOwner()[faceI];
         label bFaceI = faceI - mesh_.nInternalFaces();
 
-        if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
+        const vectorList& ownNormals = cellSurfNormals[own];
+        const labelList& ownLevels = cellSurfLevels[own];
+        const vectorList& neiNormals = neiSurfaceNormals[bFaceI];
+
+        // Special case: owner normals set is a subset of the neighbour
+        // normals. Do not do curvature refinement since other cell's normals
+        // do not add any information. Avoids weird corner extensions of
+        // refinement regions.
+        bool ownIsSubset =
+            countMatches(ownNormals, neiNormals)
+         == ownNormals.size();
+
+        bool neiIsSubset =
+            countMatches(neiNormals, ownNormals)
+         == neiNormals.size();
+
+
+        if (!ownIsSubset && !neiIsSubset)
         {
-            // Have valid data on both sides. Check curvature.
-            if ((cellMaxNormal[own] & neiBndMaxNormal[bFaceI]) < curvature)
+            // n^2 comparison of between ownNormals and neiNormals
+            for (label i = 0; !reachedLimit &&  i < ownNormals.size(); i++)
             {
-                if
-                (
-                    !markForRefine
-                    (
-                        cellMaxLevel[own],
-                        nAllowRefine,
-                        refineCell[own],
-                        nRefine
-                    )
-                )
+                for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
                 {
-                    if (debug)
+                    // Have valid data on both sides. Check curvature.
+                    if ((ownNormals[i] & neiNormals[j]) < curvature)
                     {
-                        Pout<< "Stopped refining since reaching my cell"
-                            << " limit of " << mesh_.nCells()+7*nRefine
-                            << endl;
+                        if (cellLevel[own] < ownLevels[i])
+                        {
+                            if
+                            (
+                                !markForRefine
+                                (
+                                    ownLevels[i],
+                                    nAllowRefine,
+                                    refineCell[own],
+                                    nRefine
+                                )
+                            )
+                            {
+                                if (debug)
+                                {
+                                    Pout<< "Stopped refining since reaching"
+                                        << " my cell limit of "
+                                        << mesh_.nCells()+7*nRefine
+                                        << endl;
+                                }
+                                reachedLimit = true;
+                                break;
+                            }
+                        }
                     }
-                    break;
                 }
             }
         }
     }
 
+
     if
     (
         returnReduce(nRefine, sumOp<label>())
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
index 70d9cb3bb10007ae9575aa1afd2b54f3e4148267..83e3389a5cdab446efd255076c14773dba2230dd 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
@@ -241,7 +241,7 @@ void meshRefinement::collectAndPrint
         Pstream::blocking
     );
 
-    Field<T> allData;
+    List<T> allData;
     globalPoints.gather
     (
         Pstream::worldComm,
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index bee6ada5a82a09939712052a3a318ea9fcc05f20..8d4180d6b7b7c200ca25ad9e5c8c42b042bb0410 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -720,7 +720,7 @@ void Foam::refinementSurfaces::findAllHigherIntersections
     labelList pRegions;
     vectorField pNormals;
 
-    forAll(surfaces(), surfI)
+    forAll(surfaces_, surfI)
     {
         const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
 
diff --git a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
index bb09ecb0849605e4c971d25a596738b8323d6e3c..d4699870a47a3a86f1c3ccd99988ffe4a049597e 100644
--- a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
+++ b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.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
@@ -56,12 +56,14 @@ Foam::setsToFaceZone::setsToFaceZone
 (
     const polyMesh& mesh,
     const word& faceSetName,
-    const word& cellSetName
+    const word& cellSetName,
+    const Switch& flip
 )
 :
     topoSetSource(mesh),
     faceSetName_(faceSetName),
-    cellSetName_(cellSetName)
+    cellSetName_(cellSetName),
+    flip_(flip)
 {}
 
 
@@ -74,7 +76,8 @@ Foam::setsToFaceZone::setsToFaceZone
 :
     topoSetSource(mesh),
     faceSetName_(dict.lookup("faceSet")),
-    cellSetName_(dict.lookup("cellSet"))
+    cellSetName_(dict.lookup("cellSet")),
+    flip_(dict.lookupOrDefault("flip", false))
 {}
 
 
@@ -87,7 +90,8 @@ Foam::setsToFaceZone::setsToFaceZone
 :
     topoSetSource(mesh),
     faceSetName_(checkIs(is)),
-    cellSetName_(checkIs(is))
+    cellSetName_(checkIs(is)),
+    flip_(false)
 {}
 
 
@@ -136,7 +140,7 @@ void Foam::setsToFaceZone::applyToSet
 
                 if (!fzSet.found(faceI))
                 {
-                    bool flip = false;
+                    bool flipFace = false;
 
                     label own = mesh_.faceOwner()[faceI];
                     bool ownFound = cSet.found(own);
@@ -148,11 +152,11 @@ void Foam::setsToFaceZone::applyToSet
 
                         if (ownFound && !neiFound)
                         {
-                            flip = false;
+                            flipFace = false;
                         }
                         else if (!ownFound && neiFound)
                         {
-                            flip = true;
+                            flipFace = true;
                         }
                         else
                         {
@@ -174,11 +178,17 @@ void Foam::setsToFaceZone::applyToSet
                     }
                     else
                     {
-                        flip = !ownFound;
+                        flipFace = !ownFound;
+                    }
+
+
+                    if (flip_)
+                    {
+                        flipFace = !flipFace;
                     }
 
                     newAddressing.append(faceI);
-                    newFlipMap.append(flip);
+                    newFlipMap.append(flipFace);
                 }
             }
 
diff --git a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.H b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.H
index 686849cc5000fedf3f3e3a1520272e14f41d036c..72d2f6501e1d8656b135553316ebb4e653be4275 100644
--- a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.H
+++ b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.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
@@ -36,6 +36,7 @@ SourceFiles
 #define setsToFaceZone_H
 
 #include "topoSetSource.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,10 +57,13 @@ class setsToFaceZone
         static addToUsageTable usage_;
 
         //- Name of set to use
-        word faceSetName_;
+        const word faceSetName_;
 
         //- Name of set to use
-        word cellSetName_;
+        const word cellSetName_;
+
+        //- Whether cellSet is slave cells or master cells
+        const Switch flip_;
 
 public:
 
@@ -73,7 +77,8 @@ public:
         (
             const polyMesh& mesh,
             const word& faceSetName,
-            const word& cellSetName
+            const word& cellSetName,
+            const Switch& flip
         );
 
         //- Construct from dictionary
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
index 9fd2eb7777741f08f321aede555ed40ddb7ef753..a23933030b1a519945789536653ebead7e836000 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
@@ -66,7 +66,11 @@ void thixotropicViscosity::updateMu()
     const volScalarField filmMass("filmMass", film.netMass() + mSMALL);
 
     // weighting field to blend new and existing mass contributions
-    const volScalarField w("w", max(0.0, min(1.0, deltaMass/filmMass)));
+    const volScalarField w
+    (
+        "w",
+        max(scalar(0.0), min(scalar(1.0), deltaMass/filmMass))
+    );
 
     // evaluate thixotropic viscosity
     volScalarField muThx("muThx", muInf_/(sqr(1.0 - K_*lambda_) + ROOTVSMALL));
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSchemes
index f55c52ee564855bde82231ad16e98d275780af37..0c84553dbc46a185e94d5a3566429c14cd6a9c71 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
 
     "div\(alphaPhi.*,U.*\)"     Gauss limitedLinearV 1;
     "div\(phi.*,U.*\)"          Gauss limitedLinearV 1;
-    "div\(phi.*,.*rho.*\)"      Gauss linear;
+    "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
 
     "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
     "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes
index d3f22a16e83654868ecdc438167a2bd61f8df600..0ccc0c9b67888ba13c9ec849500c9b1e08d8fc52 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
 
     "div\(alphaPhi.*,U.*\)"     Gauss limitedLinearV 1;
     "div\(phi.*,U.*\)"          Gauss limitedLinearV 1;
-    "div\(phi.*,.*rho.*\)"      Gauss linear;
+    "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
 
     "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
     "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSchemes
index 2546db411275847e1db31badc51db945562292ac..1c6548cac37c3d48da4c96941dfb5f51b1683442 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
 
     "div\(alphaPhi.*,U.*\)"     Gauss limitedLinearV 1;
     "div\(phi.*,U.*\)"          Gauss limitedLinearV 1;
-    "div\(phi.*,.*rho.*\)"      Gauss linear;
+    "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
 
     "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
     "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSchemes
index 549178c4c6d23cc27d375e5f3ec2f1c7084672ee..0bd5b9026ec94c19b121fe86056002cf740600e8 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
 
     "div\(alphaPhi.*,U.*\)"     Gauss limitedLinearV 1;
     "div\(phi.*,U.*\)"          Gauss limitedLinearV 1;
-    "div\(phi.*,.*rho.*\)"      Gauss linear;
+    "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
 
     "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
     "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution
index 6e38c9cc85f6257bff5c7408a589a29ee25271f3..a25ccff8c62a0d6c8289dc37ad95385ffaaaf771 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution
@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    alpha
+    alpha.air
     {
         nAlphaCorr      1;
         nAlphaSubCycles 2;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes
index 3fd8876cb600ec235534eb44522c437e989fb8df..c405673fdaadf30aa5668d1f433553970e5ebf59 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes
@@ -36,7 +36,7 @@ divSchemes
 
     "div\(alphaPhi.*,U.*\)"     Gauss limitedLinearV 1;
     "div\(phi.*,U.*\)"          Gauss limitedLinearV 1;
-    "div\(phi.*,.*rho.*\)"      Gauss linear;
+    "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
 
     "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
     "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSchemes
index 9fc65bf8cbc74f0036a8175ef25d5f528fb57565..0360d1b8ce1d847bad66314a1e19d67cd5add7c7 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
 
     "div\(alphaPhi.*,U.*\)"     Gauss limitedLinearV 1;
     "div\(phi.*,U.*\)"          Gauss limitedLinearV 1;
-    "div\(phi.*,.*rho.*\)"      Gauss linear;
+    "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
 
     "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
     "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSchemes
index c29b10072aedbfddf0ac04218b035914f8091684..61f5962400d0546881c35db9f00fdd5ec799191d 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
 
     "div\(alphaPhi.*,U.*\)"     Gauss limitedLinearV 1;
     "div\(phi.*,U.*\)"          Gauss limitedLinearV 1;
-    "div\(phi.*,.*rho.*\)"      Gauss linear;
+    "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
 
     "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
     "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;