diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 9b09d73bbe..858d83749e 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    rhoSimpleFoam
+    rhoPimpleFoam
 
 Description
     Transient solver for turbulent flow of compressible fluids for
diff --git a/applications/solvers/incompressible/pimpleFoam/Make/files b/applications/solvers/incompressible/pimpleFoam/Make/files
index 92943d3370..f4ec4391df 100644
--- a/applications/solvers/incompressible/pimpleFoam/Make/files
+++ b/applications/solvers/incompressible/pimpleFoam/Make/files
@@ -1,3 +1,3 @@
 pimpleFoam.C
 
-EXE = $(FOAM_USER_APPBIN)/pimpleFoam
+EXE = $(FOAM_APPBIN)/pimpleFoam
diff --git a/applications/test/speed/Make/options b/applications/test/speed/Make/options
index e69de29bb2..e680561987 100644
--- a/applications/test/speed/Make/options
+++ b/applications/test/speed/Make/options
@@ -0,0 +1 @@
+EXE_INC = /* -ffast-math -mtune=core2 */
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 1df11e06da..4737f38e53 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -370,7 +370,7 @@ int main(int argc, char *argv[])
         // Refinement parameters
         refinementParameters refineParams(refineDict);
 
-        refineDriver.doRefine(refineDict, refineParams, wantSnap);
+        refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict);
 
         writeMesh
         (
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
index 48f52d7880..c07b950d0a 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
@@ -491,6 +491,8 @@ void Foam::autoHexMeshDriver::doMesh()
 
     if (wantRefine)
     {
+        const dictionary& motionDict = dict_.subDict("motionDict");
+
         autoRefineDriver refineDriver
         (
             meshRefinerPtr_(),
@@ -502,7 +504,7 @@ void Foam::autoHexMeshDriver::doMesh()
         // Get all the refinement specific params
         refinementParameters refineParams(dict_, 1);
 
-        refineDriver.doRefine(dict_, refineParams, wantSnap);
+        refineDriver.doRefine(dict_, refineParams, wantSnap, motionDict);
 
         // Write mesh
         writeMesh("Refined mesh");
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index 7bbe58ff8f..ca0c0cbddc 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -497,7 +497,8 @@ Foam::label Foam::autoRefineDriver::shellRefine
 void Foam::autoRefineDriver::baffleAndSplitMesh
 (
     const refinementParameters& refineParams,
-    const bool handleSnapProblems
+    const bool handleSnapProblems,
+    const dictionary& motionDict
 )
 {
     Info<< nl
@@ -514,6 +515,7 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
     (
         handleSnapProblems,
         !handleSnapProblems,            // merge free standing baffles?
+        motionDict,
         const_cast<Time&>(mesh.time()),
         globalToPatch_,
         refineParams.keepPoints()[0]
@@ -568,7 +570,8 @@ void Foam::autoRefineDriver::zonify
 void Foam::autoRefineDriver::splitAndMergeBaffles
 (
     const refinementParameters& refineParams,
-    const bool handleSnapProblems
+    const bool handleSnapProblems,
+    const dictionary& motionDict
 )
 {
     Info<< nl
@@ -588,6 +591,7 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
     (
         handleSnapProblems,
         false,                  // merge free standing baffles?
+        motionDict,
         const_cast<Time&>(mesh.time()),
         globalToPatch_,
         refineParams.keepPoints()[0]
@@ -685,7 +689,8 @@ void Foam::autoRefineDriver::doRefine
 (
     const dictionary& refineDict,
     const refinementParameters& refineParams,
-    const bool prepareForSnapping
+    const bool prepareForSnapping,
+    const dictionary& motionDict
 )
 {
     Info<< nl
@@ -734,13 +739,13 @@ void Foam::autoRefineDriver::doRefine
     );
 
     // Introduce baffles at surface intersections
-    baffleAndSplitMesh(refineParams, prepareForSnapping);
+    baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
 
     // Mesh is at its finest. Do optional zoning.
     zonify(refineParams);
 
     // Pull baffles apart
-    splitAndMergeBaffles(refineParams, prepareForSnapping);
+    splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
 
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
index 8aa33855de..f1747a29be 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
@@ -46,7 +46,6 @@ namespace Foam
 class featureEdgeMesh;
 class refinementParameters;
 
-
 /*---------------------------------------------------------------------------*\
                            Class autoRefineDriver Declaration
 \*---------------------------------------------------------------------------*/
@@ -112,7 +111,8 @@ class autoRefineDriver
         void baffleAndSplitMesh
         (
             const refinementParameters& refineParams,
-            const bool handleSnapProblems
+            const bool handleSnapProblems,
+            const dictionary& motionDict
         );
 
         //- Add zones
@@ -121,7 +121,8 @@ class autoRefineDriver
         void splitAndMergeBaffles
         (
             const refinementParameters& refineParams,
-            const bool handleSnapProblems
+            const bool handleSnapProblems,
+            const dictionary& motionDict
         );
 
         //- Merge refined boundary faces (from exposing coarser cell)
@@ -163,7 +164,8 @@ public:
         (
             const dictionary& refineDict,
             const refinementParameters& refineParams,
-            const bool prepareForSnapping
+            const bool prepareForSnapping,
+            const dictionary& motionDict
         );
 };
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index cc866d13ca..520e6f69ef 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -1294,6 +1294,169 @@ void Foam::autoSnapDriver::scaleMesh
 }
 
 
+// After snapping: correct patching according to nearest surface.
+// Code is very similar to calcNearestSurface.
+// - calculate face-wise snap distance as max of point-wise
+// - calculate face-wise nearest surface point
+// - repatch face according to patch for surface point.
+Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
+(
+    const snapParameters& snapParams,
+    const labelList& adaptPatchIDs
+)
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
+
+    Info<< "Repatching faces according to nearest surface ..." << endl;
+
+    // Get the labels of added patches.
+    autoPtr<indirectPrimitivePatch> ppPtr
+    (
+        meshRefinement::makePatch
+        (
+            mesh,
+            adaptPatchIDs
+        )
+    );
+    indirectPrimitivePatch& pp = ppPtr();
+
+    // Divide surfaces into zoned and unzoned
+    labelList zonedSurfaces;
+    labelList unzonedSurfaces;
+    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
+
+
+    // Faces that do not move
+    PackedList<1> isZonedFace(mesh.nFaces(), 0);
+    {
+        // 1. All faces on zoned surfaces
+        const wordList& faceZoneNames = surfaces.faceZoneNames();
+        const faceZoneMesh& fZones = mesh.faceZones();
+
+        forAll(zonedSurfaces, i)
+        {
+            label zoneSurfI = zonedSurfaces[i];
+
+            label zoneI = fZones.findZoneID(faceZoneNames[zoneSurfI]);
+
+            const faceZone& fZone = fZones[zoneI];
+
+            forAll(fZone, i)
+            {
+                isZonedFace.set(fZone[i], 1);
+            }  
+        }
+    }
+
+
+    // Determine per pp face which patch it should be in
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Patch that face should be in
+    labelList closestPatch(pp.size(), -1);
+    {
+        // face snap distance as max of point snap distance
+        scalarField faceSnapDist(pp.size(), -GREAT);
+        {
+            // Distance to attract to nearest feature on surface
+            const scalarField snapDist(calcSnapDistance(snapParams, pp));
+
+            const faceList& localFaces = pp.localFaces();
+
+            forAll(localFaces, faceI)
+            {
+                const face& f = localFaces[faceI];
+
+                forAll(f, fp)
+                {
+                    faceSnapDist[faceI] = max
+                    (
+                        faceSnapDist[faceI],
+                        snapDist[f[fp]]
+                    );
+                }
+            }
+        }
+
+        pointField localFaceCentres(pp.size());
+        forAll(pp, i)
+        {
+            localFaceCentres[i] = mesh.faceCentres()[pp.addressing()[i]];
+        }
+
+        // Get nearest surface and region
+        labelList hitSurface;
+        labelList hitRegion;
+        surfaces.findNearestRegion
+        (
+            unzonedSurfaces,
+            localFaceCentres,
+            sqr(4*faceSnapDist),    // sqr of attract distance
+            hitSurface,
+            hitRegion
+        );
+
+        // Get patch
+        forAll(pp, i)
+        {
+            label faceI = pp.addressing()[i];
+
+            if (hitSurface[i] != -1 && (isZonedFace.get(faceI) == 0))
+            {
+                closestPatch[i] = globalToPatch_
+                [
+                    surfaces.globalRegion
+                    (
+                        hitSurface[i],
+                        hitRegion[i]
+                    )
+                ];
+            }
+        }
+    }
+
+
+    // Change those faces for which there is a different closest patch
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList ownPatch(mesh.nFaces(), -1);
+    labelList neiPatch(mesh.nFaces(), -1);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        forAll(pp, i)
+        {
+            ownPatch[pp.start()+i] = patchI;
+            neiPatch[pp.start()+i] = patchI;
+        }
+    }
+
+    label nChanged = 0;
+    forAll(closestPatch, i)
+    {
+        label faceI = pp.addressing()[i];
+
+        if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[faceI])
+        {
+            ownPatch[faceI] = closestPatch[i];
+            neiPatch[faceI] = closestPatch[i];
+            nChanged++;
+        }
+    }
+
+    Info<< "Repatched " << returnReduce(nChanged, sumOp<label>())
+        << " faces in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl
+        << endl;
+
+    return meshRefiner_.createBaffles(ownPatch, neiPatch);
+}
+
+
 void Foam::autoSnapDriver::doSnap
 (
     const dictionary& snapDict,
@@ -1329,7 +1492,7 @@ void Foam::autoSnapDriver::doSnap
         );
         indirectPrimitivePatch& pp = ppPtr();
 
-        // Distance to attact to nearest feature on surface
+        // Distance to attract to nearest feature on surface
         const scalarField snapDist(calcSnapDistance(snapParams, pp));
 
 
@@ -1383,6 +1546,9 @@ void Foam::autoSnapDriver::doSnap
 
     // Merge any introduced baffles.
     mergeZoneBaffles(baffles);
+
+    // Repatch faces according to nearest.
+    repatchToSurface(snapParams, adaptPatchIDs);
 }
 
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 82cf448b98..1604e35859 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -215,6 +215,14 @@ public:
                 motionSmoother&
             );
 
+            //- Repatch faces according to surface nearest the face centre
+            autoPtr<mapPolyMesh> repatchToSurface
+            (
+                const snapParameters& snapParams,
+                const labelList& adaptPatchIDs
+            );
+
+
             void doSnap
             (
                 const dictionary& snapDict,
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 4fe3005766..7f5bc29a79 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -50,6 +50,7 @@ License
 #include "globalIndex.H"
 #include "meshTools.H"
 #include "OFstream.H"
+#include "geomDecomp.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -426,15 +427,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
 
 // Determine for multi-processor regions the lowest numbered cell on the lowest
 // numbered processor.
-void Foam::meshRefinement::getRegionMaster
+void Foam::meshRefinement::getCoupledRegionMaster
 (
+    const globalIndex& globalCells,
     const boolList& blockedFace,
     const regionSplit& globalRegion,
     Map<label>& regionToMaster
 ) const
 {
-    globalIndex globalCells(mesh_.nCells());
-
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
     forAll(patches, patchI)
@@ -483,6 +483,274 @@ void Foam::meshRefinement::getRegionMaster
 }
 
 
+void Foam::meshRefinement::calcLocalRegions
+(
+    const globalIndex& globalCells,
+    const labelList& globalRegion,
+    const Map<label>& coupledRegionToMaster,
+
+    Map<label>& globalToLocalRegion,
+    pointField& localPoints
+) const
+{
+    globalToLocalRegion.resize(globalRegion.size());
+    DynamicList<point> localCc(globalRegion.size()/2);
+
+    forAll(globalRegion, cellI)
+    {
+        Map<label>::const_iterator fndMaster =
+            coupledRegionToMaster.find(globalRegion[cellI]);
+
+        if (fndMaster != coupledRegionToMaster.end())
+        {
+            // Multi-processor region.
+            if (globalCells.toGlobal(cellI) == fndMaster())
+            {
+                // I am master. Allocate region for me.
+                globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
+                localCc.append(mesh_.cellCentres()[cellI]);
+            }
+        }
+        else
+        {
+            // Single processor region.
+            if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
+            {
+                localCc.append(mesh_.cellCentres()[cellI]);
+            }
+        }
+    }
+    localCc.shrink();
+    localPoints.transfer(localCc);
+
+    if (localPoints.size() != globalToLocalRegion.size())
+    {
+        FatalErrorIn("calcLocalRegions(..)")
+            << "localPoints:" << localPoints.size()
+            << " globalToLocalRegion:" << globalToLocalRegion.size()
+            << exit(FatalError);
+    }
+}
+
+
+Foam::label Foam::meshRefinement::getShiftedRegion
+(
+    const globalIndex& indexer,
+    const Map<label>& globalToLocalRegion,
+    const Map<label>& coupledRegionToShifted,
+    const label globalRegion
+)
+{
+    Map<label>::const_iterator iter =
+        globalToLocalRegion.find(globalRegion);
+
+    if (iter != globalToLocalRegion.end())
+    {
+        // Region is 'owned' locally. Convert local region index into global.
+        return indexer.toGlobal(iter());
+    }
+    else
+    {
+        return coupledRegionToShifted[globalRegion];
+    }
+}
+
+
+// Add if not yet present
+void Foam::meshRefinement::addUnique(const label elem, labelList& lst)
+{
+    if (findIndex(lst, elem) == -1)
+    {
+        label sz = lst.size();
+        lst.setSize(sz+1);
+        lst[sz] = elem;
+    }
+}
+
+
+void Foam::meshRefinement::calcRegionRegions
+(
+    const labelList& globalRegion,
+    const Map<label>& globalToLocalRegion,
+    const Map<label>& coupledRegionToMaster,
+    labelListList& regionRegions
+) const
+{
+    // Global region indexing since we now know the shifted regions.
+    globalIndex shiftIndexer(globalToLocalRegion.size());
+
+    // Redo the coupledRegionToMaster to be in shifted region indexing.
+    Map<label> coupledRegionToShifted(coupledRegionToMaster.size());
+    forAllConstIter(Map<label>, coupledRegionToMaster, iter)
+    {
+        label region = iter.key();
+
+        Map<label>::const_iterator fndRegion = globalToLocalRegion.find(region);
+
+        if (fndRegion != globalToLocalRegion.end())
+        {
+            // A local cell is master of this region. Get its shifted region.
+            coupledRegionToShifted.insert
+            (
+                region,
+                shiftIndexer.toGlobal(fndRegion())
+            );
+        }
+        Pstream::mapCombineGather(coupledRegionToShifted, minEqOp<label>());
+        Pstream::mapCombineScatter(coupledRegionToShifted);
+    }
+
+
+    // Determine region-region connectivity.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // This is for all my regions (so my local ones or the ones I am master
+    // of) the neighbouring region indices.
+
+
+    // Transfer lists.
+    PtrList<HashSet<edge, Hash<edge> > > regionConnectivity(Pstream::nProcs());
+    forAll(regionConnectivity, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            regionConnectivity.set
+            (
+                procI,
+                new HashSet<edge, Hash<edge> >
+                (
+                    coupledRegionToShifted.size()
+                  / Pstream::nProcs()
+                )
+            );
+        }
+    }
+
+
+    // Connectivity. For all my local regions the connected regions.
+    regionRegions.setSize(globalToLocalRegion.size());
+
+    // Add all local connectivity to regionRegions; add all non-local
+    // to the transferlists.
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        label ownRegion = globalRegion[mesh_.faceOwner()[faceI]];
+        label neiRegion = globalRegion[mesh_.faceNeighbour()[faceI]];
+
+        if (ownRegion != neiRegion)
+        {
+            label shiftOwnRegion = getShiftedRegion
+            (
+                shiftIndexer,
+                globalToLocalRegion,
+                coupledRegionToShifted,
+                ownRegion
+            );
+            label shiftNeiRegion = getShiftedRegion
+            (
+                shiftIndexer,
+                globalToLocalRegion,
+                coupledRegionToShifted,
+                neiRegion
+            );
+
+
+            // Connection between two regions. Send to owner of region.
+            // - is ownRegion 'owned' by me
+            // - is neiRegion 'owned' by me
+
+            if (shiftIndexer.isLocal(shiftOwnRegion))
+            {
+                label localI = shiftIndexer.toLocal(shiftOwnRegion);
+                addUnique(shiftNeiRegion, regionRegions[localI]);
+            }
+            else
+            {
+                label masterProc = shiftIndexer.whichProcID(shiftOwnRegion);
+                regionConnectivity[masterProc].insert
+                (
+                    edge(shiftOwnRegion, shiftNeiRegion)
+                );
+            }
+
+            if (shiftIndexer.isLocal(shiftNeiRegion))
+            {
+                label localI = shiftIndexer.toLocal(shiftNeiRegion);
+                addUnique(shiftOwnRegion, regionRegions[localI]);
+            }
+            else
+            {
+                label masterProc = shiftIndexer.whichProcID(shiftNeiRegion);
+                regionConnectivity[masterProc].insert
+                (
+                    edge(shiftOwnRegion, shiftNeiRegion)
+                );
+            }
+        }
+    }
+
+
+    // Send
+    forAll(regionConnectivity, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            OPstream str(Pstream::blocking, procI);
+            str << regionConnectivity[procI];
+        }
+    }
+    // Receive
+    forAll(regionConnectivity, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            IPstream str(Pstream::blocking, procI);
+            str >> regionConnectivity[procI];
+        }
+    }
+
+    // Add to addressing.
+    forAll(regionConnectivity, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            for
+            (
+                HashSet<edge, Hash<edge> >::const_iterator iter =
+                    regionConnectivity[procI].begin();
+                iter != regionConnectivity[procI].end();
+                ++iter
+            )
+            {
+                const edge& e = iter.key();
+
+                bool someLocal = false;
+                if (shiftIndexer.isLocal(e[0]))
+                {
+                    label localI = shiftIndexer.toLocal(e[0]);
+                    addUnique(e[1], regionRegions[localI]);
+                    someLocal = true;
+                }
+                if (shiftIndexer.isLocal(e[1]))
+                {
+                    label localI = shiftIndexer.toLocal(e[1]);
+                    addUnique(e[0], regionRegions[localI]);
+                    someLocal = true;
+                }
+
+                if (!someLocal)
+                {
+                    FatalErrorIn("calcRegionRegions(..)")
+                        << "Received from processor " << procI
+                        << " connection " << e
+                        << " where none of the elements is local to me."
+                        << abort(FatalError);
+                }
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
@@ -598,88 +866,88 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
     // the region might span multiple domains so we want to get
     // a "region master" per domain. Note that multi-processor
     // regions can only occur on cells on coupled patches.
-
-
-    // Determine per coupled region the master cell (lowest numbered cell
-    // on lowest numbered processor)
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    Map<label> regionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
-    getRegionMaster(blockedFace, globalRegion, regionToMaster);
+    // Note: since the number of regions does not change by this the
+    // process can be seen as just a shift of a region onto a single
+    // processor.
 
 
     // Global cell numbering engine
     globalIndex globalCells(mesh_.nCells());
 
 
-    // Determine cell centre per region
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Now we can divide regions into
-    // - single-processor regions (almost all)
-    //   - allocate local region + coordinate for it
-    // - multi-processor for which I am master
-    //   - allocate local region + coordinate for it
-    // - multi-processor for which I am not master
-    //   - do not allocate region.
-    //   - but inherit new distribution from master.
+    // Determine per coupled region the master cell (lowest numbered cell
+    // on lowest numbered processor)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // (does not determine master for non-coupled=fully-local regions)
 
-    Map<label> globalToLocalRegion(mesh_.nCells());
-    DynamicList<point> localCc(mesh_.nCells()/10);
+    Map<label> coupledRegionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
+    getCoupledRegionMaster
+    (
+        globalCells,
+        blockedFace,
+        globalRegion,
+        coupledRegionToMaster
+    );
 
-    forAll(globalRegion, cellI)
-    {
-        Map<label>::const_iterator fndMaster =
-            regionToMaster.find(globalRegion[cellI]);
+    // Determine my regions
+    // ~~~~~~~~~~~~~~~~~~~~
+    // These are the fully local ones or the coupled ones of which I am master.
 
-        if (fndMaster != regionToMaster.end())
-        {
-            // Multi-processor region.
-            if (globalCells.toGlobal(cellI) == fndMaster())
-            {
-                // I am master. Allocate region for me.
-                globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
-                localCc.append(mesh_.cellCentres()[cellI]);
-            }
-        }
-        else
-        {
-            // Single processor region.
-            Map<label>::iterator iter =
-                globalToLocalRegion.find(globalRegion[cellI]);
-
-            if (iter == globalToLocalRegion.end())
-            {
-                globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
-                localCc.append(mesh_.cellCentres()[cellI]);
-            }
-        }
-    }
-    localCc.shrink();
+    Map<label> globalToLocalRegion;
     pointField localPoints;
-    localPoints.transfer(localCc);
+    calcLocalRegions
+    (
+        globalCells,
+        globalRegion,
+        coupledRegionToMaster,
+
+        globalToLocalRegion,
+        localPoints
+    );
 
 
-    // Call decomposer on localCc
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    labelList localDistribution = decomposer.decompose(localPoints);
+    // Find distribution for regions
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+    labelList regionDistribution;
 
-    // Distribute the destination processor for coupled regions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    Map<label> regionToDist(regionToMaster.size());
-
-    forAllConstIter(Map<label>, regionToMaster, iter)
+    if (isA<geomDecomp>(decomposer))
     {
-        if (globalCells.isLocal(iter()))
+        regionDistribution = decomposer.decompose(localPoints);
+    }
+    else
+    {
+        labelListList regionRegions;
+        calcRegionRegions
+        (
+            globalRegion,
+            globalToLocalRegion,
+            coupledRegionToMaster,
+
+            regionRegions
+        );
+
+        regionDistribution = decomposer.decompose(regionRegions, localPoints);
+    }
+
+
+
+    // Convert region-based decomposition back to cell-based one
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Transfer destination processor back to all. Use global reduce for now.
+    Map<label> regionToDist(coupledRegionToMaster.size());
+    forAllConstIter(Map<label>, coupledRegionToMaster, iter)
+    {
+        label region = iter.key();
+
+        Map<label>::const_iterator regionFnd = globalToLocalRegion.find(region);
+
+        if (regionFnd != globalToLocalRegion.end())
         {
-            // Master cell is local.
-            regionToDist.insert
-            (
-                iter.key(),
-                localDistribution[globalToLocalRegion[iter.key()]]
-            );
+            // Master cell is local. Store my distribution.
+            regionToDist.insert(iter.key(), regionDistribution[regionFnd()]);
         }
         else
         {
@@ -691,29 +959,27 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
     Pstream::mapCombineScatter(regionToDist);
 
 
-    // Convert region-based decomposition back to cell-based one
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+    // Determine destination for all cells
     labelList distribution(mesh_.nCells());
 
     forAll(globalRegion, cellI)
     {
-        Map<label>::const_iterator fndMaster =
+        Map<label>::const_iterator fndRegion =
             regionToDist.find(globalRegion[cellI]);
 
-        if (fndMaster != regionToDist.end())
+        if (fndRegion != regionToDist.end())
         {
-            // Special handling
-            distribution[cellI] = fndMaster();
+            distribution[cellI] = fndRegion();
         }
         else
         {
             // region is local to the processor.
             label localRegionI = globalToLocalRegion[globalRegion[cellI]];
 
-            distribution[cellI] = localDistribution[localRegionI];
+            distribution[cellI] = regionDistribution[localRegionI];
         }
     }
+
     return distribution;
 }
 
@@ -839,9 +1105,18 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
         if (debug)
         {
-            Pout<< "Wanted distribution:"
-                << distributor.countCells(distribution)
-                << endl;
+            labelList nProcCells(distributor.countCells(distribution));
+            Pout<< "Wanted distribution:" << nProcCells << endl;
+
+            Pstream::listCombineGather(nProcCells, plusEqOp<label>());
+            Pstream::listCombineScatter(nProcCells);
+
+            Pout<< "Wanted resulting decomposition:" << endl;
+            forAll(nProcCells, procI)
+            {
+                Pout<< "    " << procI << '\t' << nProcCells[procI] << endl;
+            }
+            Pout<< endl;
         }
         // Do actual sending/receiving of mesh
         map = distributor.distribute(distribution);
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 012ef08d00..d3eb52b7c4 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -68,6 +68,7 @@ class featureEdgeMesh;
 class fvMeshDistribute;
 class searchableSurface;
 class regionSplit;
+class globalIndex;
 
 /*---------------------------------------------------------------------------*\
                            Class meshRefinement Declaration
@@ -157,15 +158,6 @@ private:
             labelList& neiLevel,
             pointField& neiCc
         ) const;
-        ////- Calculate coupled boundary end vector and refinement level. Sort
-        ////  so both sides of coupled face have data in same order.
-        //void calcCanonicalBoundaryData
-        //(
-        //    labelList& leftLevel,
-        //    pointField& leftCc,
-        //    labelList& rightLevel,
-        //    pointField& rightCc
-        //)  const;
 
         //- Find any intersection of surface. Store in surfaceIndex_.
         void updateIntersections(const labelList& changedFaces);
@@ -194,15 +186,51 @@ private:
             const label exposedPatchI
         );
 
-        //- Used in decomposeCombineRegions. Given global region per cell
-        //  determines master processor/cell for regions straddling
-        //  procboundaries.
-        void getRegionMaster
-        (
-            const boolList& blockedFace,
-            const regionSplit& globalRegion,
-            Map<label>& regionToMaster
-        ) const;
+        // For decomposeCombineRegions
+
+            //- Used in decomposeCombineRegions. Given global region per cell
+            //  determines master processor/cell for regions straddling
+            //  procboundaries.
+            void getCoupledRegionMaster
+            (
+                const globalIndex& globalCells,
+                const boolList& blockedFace,
+                const regionSplit& globalRegion,
+                Map<label>& regionToMaster
+            ) const;
+
+            //- Determine regions that are local to me or coupled ones that
+            //  are owned by me. Determine representative location.
+            void calcLocalRegions
+            (
+                const globalIndex& globalCells,
+                const labelList& globalRegion,
+                const Map<label>& coupledRegionToMaster,
+
+                Map<label>& globalToLocalRegion,
+                pointField& localPoints
+            ) const;
+
+            //- Convert region into global index.
+            static label getShiftedRegion
+            (
+                const globalIndex& indexer,
+                const Map<label>& globalToLocalRegion,
+                const Map<label>& coupledRegionToShifted,
+                const label globalRegion
+            );
+
+            //- helper: add element if not in list. Linear search.
+            static void addUnique(const label, labelList&);
+
+            //- Calculate region connectivity. Major communication.
+            void calcRegionRegions
+            (
+                const labelList& globalRegion,
+                const Map<label>& globalToLocalRegion,
+                const Map<label>& coupledRegionToMaster,
+                labelListList& regionRegions
+            ) const;
 
 
         // Refinement candidate selection
@@ -334,6 +362,14 @@ private:
                 const labelList& globalToPatch
             ) const;
 
+            //- Returns list with for every internal face -1 or the patch
+            //  they should be baffled into.
+            labelList markFacesOnProblemCellsGeometric
+            (
+                const dictionary& motionDict,
+                const labelList& globalToPatch
+            ) const;
+
 
         // Baffle merging
 
@@ -552,6 +588,7 @@ public:
             (
                 const bool handleSnapProblems,
                 const bool mergeFreeStanding,
+                const dictionary& motionDict,
                 Time& runTime,
                 const labelList& globalToPatch,
                 const point& keepPoint
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 8b0236a7e2..e6bd37f670 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -44,9 +44,9 @@ License
 #include "OFstream.H"
 #include "regionSplit.H"
 #include "removeCells.H"
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+#include "motionSmoother.H"
+#include "polyMeshGeometry.H"
+#include "IOmanip.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -555,28 +555,28 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     boolList isBoundaryFace(mesh_.nFaces(), false);
 
     // Fill boundary data. All elements on meshed patches get marked.
-    forAll(globalToPatch, i)
+    // Get the labels of added patches.
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    forAll(adaptPatchIDs, i)
     {
-        label patchI = globalToPatch[i];
+        label patchI = adaptPatchIDs[i];
 
-        if (patchI != -1)
+        const polyPatch& pp = patches[patchI];
+
+        label faceI = pp.start();
+
+        forAll(pp, j)
         {
-            const polyPatch& pp = patches[patchI];
+            markBoundaryFace
+            (
+                faceI,
+                isBoundaryFace,
+                isBoundaryEdge,
+                isBoundaryPoint
+            );
 
-            label faceI = pp.start();
-
-            forAll(pp, j)
-            {
-                markBoundaryFace
-                (
-                    faceI,
-                    isBoundaryFace,
-                    isBoundaryEdge,
-                    isBoundaryPoint
-                );
-
-                faceI++;
-            }
+            faceI++;
         }
     }
 
@@ -872,6 +872,264 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
     return facePatch;
 }
+//XXXXXXXXXXXXXX
+// Mark faces to be baffled to prevent snapping problems. Does
+// test to find nearest surface and checks which faces would get squashed.
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
+(
+    const dictionary& motionDict,
+    const labelList& globalToPatch
+) const
+{
+    // Get the labels of added patches.
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    // Construct addressing engine.
+    autoPtr<indirectPrimitivePatch> ppPtr
+    (
+        meshRefinement::makePatch
+        (
+            mesh_,
+            adaptPatchIDs
+        )
+    );
+    const indirectPrimitivePatch& pp = ppPtr();
+    const pointField& localPoints = pp.localPoints();
+    const labelList& meshPoints = pp.meshPoints();
+
+    // Find nearest (non-baffle) surface
+    pointField newPoints(mesh_.points());
+    {
+        List<pointIndexHit> hitInfo;
+        labelList hitSurface;
+        surfaces_.findNearest
+        (
+            surfaces_.getUnnamedSurfaces(),
+            localPoints,
+            scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
+            hitSurface,
+            hitInfo
+        );
+    
+        forAll(hitInfo, i)
+        {
+            if (hitInfo[i].hit())
+            {
+                //label pointI = meshPoints[i];
+                //Pout<< "   " << pointI << " moved from "
+                //    << mesh_.points()[pointI] << " by "
+                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
+                //    << endl;
+                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
+            }
+        }
+    }
+
+    // Per face (internal or coupled!) the patch that the
+    // baffle should get (or -1).
+    labelList facePatch(mesh_.nFaces(), -1);
+    // Count of baffled faces
+    label nBaffleFaces = 0;
+
+
+//    // Sync position? Or not since same face on both side so just sync
+//    // result of baffle.
+//
+//    const scalar minArea(readScalar(motionDict.lookup("minArea")));
+//
+//    Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
+//        << minArea << endl;
+//
+//    pointField facePoints;
+//    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+//    {
+//        const face& f = mesh_.faces()[faceI];
+//
+//        bool usesPatchPoint = false;
+//
+//        facePoints.setSize(f.size());
+//        forAll(f, fp)
+//        {
+//            Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
+//
+//            if (iter != pp.meshPointMap().end())
+//            {
+//                facePoints[fp] = newPosition[iter()];
+//                usesPatchPoint = true;
+//            }
+//            else
+//            {
+//                facePoints[fp] = mesh_.points()[f[fp]];
+//            }
+//        }
+//
+//        if (usesPatchPoint)
+//        {
+//            // Check area of face wrt original area
+//            face identFace(identity(f.size()));
+//
+//            if (identFace.mag(facePoints) < minArea)
+//            {
+//                facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                nBaffleFaces++;
+//            }
+//        }
+//    }
+//
+//
+//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+//    forAll(patches, patchI)
+//    {
+//        const polyPatch& pp = patches[patchI];
+//
+//        if (pp.coupled())
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.start()+i;
+//
+//                const face& f = mesh_.faces()[faceI];
+//
+//                bool usesPatchPoint = false;
+//
+//                facePoints.setSize(f.size());
+//                forAll(f, fp)
+//                {
+//                    Map<label>::const_iterator iter =
+//                        pp.meshPointMap().find(f[fp]);
+//
+//                    if (iter != pp.meshPointMap().end())
+//                    {
+//                        facePoints[fp] = newPosition[iter()];
+//                        usesPatchPoint = true;
+//                    }
+//                    else
+//                    {
+//                        facePoints[fp] = mesh_.points()[f[fp]];
+//                    }
+//                }
+//
+//                if (usesPatchPoint)
+//                {
+//                    // Check area of face wrt original area
+//                    face identFace(identity(f.size()));
+//
+//                    if (identFace.mag(facePoints) < minArea)
+//                    {
+//                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                        nBaffleFaces++;
+//                    }
+//                }
+//            }
+//        }
+//    }
+
+    {
+        pointField oldPoints(mesh_.points());
+        mesh_.movePoints(newPoints);
+        faceSet wrongFaces(mesh_, "wrongFaces", 100);
+        {
+            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+
+            // Just check the errors from squashing
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            const labelList allFaces(identity(mesh_.nFaces()));
+            label nWrongFaces = 0;
+
+            scalar minArea(readScalar(motionDict.lookup("minArea")));
+            if (minArea > -SMALL)
+            {
+                polyMeshGeometry::checkFaceArea
+                (
+                    false,
+                    minArea,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces with area < "
+                    << setw(5) << minArea
+                    << " m^2                            : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+
+//            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+            scalar minDet = 0.01;
+            if (minDet > -1)
+            {
+                polyMeshGeometry::checkCellDeterminant
+                (
+                    false,
+                    minDet,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    polyMeshGeometry::affectedCells(mesh_, allFaces),
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces on cells with determinant < "
+                    << setw(5) << minDet << "                : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+        }
+
+
+        forAllConstIter(faceSet, wrongFaces, iter)
+        {
+            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+
+            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+            {
+                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+                nBaffleFaces++;
+
+                //Pout<< "    " << iter.key()
+                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+                //    << " is destined for patch " << facePatch[iter.key()]
+                //    << endl;
+            }
+        }
+        // Restore points.
+        mesh_.movePoints(oldPoints);
+    }
+
+
+    Info<< "markFacesOnProblemCellsGeometric : marked "
+        << returnReduce(nBaffleFaces, sumOp<label>())
+        << " additional internal and coupled faces"
+        << " to be converted into baffles." << endl;
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        facePatch,
+        maxEqOp<label>(),
+        false               // no separation
+    );
+
+    return facePatch;
+}
+//XXXXXXXX
 
 
 // Return a list of coupled face pairs, i.e. faces that use the same vertices.
@@ -1554,6 +1812,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
 (
     const bool handleSnapProblems,
     const bool mergeFreeStanding,
+    const dictionary& motionDict,
     Time& runTime,
     const labelList& globalToPatch,
     const point& keepPoint
@@ -1584,6 +1843,12 @@ void Foam::meshRefinement::baffleAndSplitMesh
         ownPatch,
         neiPatch
     );
+
+    if (debug)
+    {
+        runTime++;
+    }
+
     createBaffles(ownPatch, neiPatch);
 
     if (debug)
@@ -1621,14 +1886,65 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
         labelList facePatch
         (
-            markFacesOnProblemCells
+            //markFacesOnProblemCells
+            //(
+            //    globalToPatch
+            //)
+            markFacesOnProblemCellsGeometric
             (
+                motionDict,
                 globalToPatch
             )
         );
         Info<< "Analyzed problem cells in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
 
+        if (debug)
+        {
+            // Dump all these faces to a faceSet.
+            faceSet problemGeom(mesh_, "problemFacesGeom", 100);
+
+            const labelList facePatchGeom
+            (
+                markFacesOnProblemCellsGeometric
+                (
+                    motionDict,
+                    globalToPatch
+                )
+            );
+            forAll(facePatchGeom, faceI)
+            {
+                if (facePatchGeom[faceI] != -1)
+                {
+                    problemGeom.insert(faceI);
+                }
+            }
+            Pout<< "Dumping " << problemGeom.size()
+                << " problem faces to " << problemGeom.objectPath() << endl;
+            problemGeom.write();
+
+
+            faceSet problemTopo(mesh_, "problemFacesTopo", 100);
+
+            const labelList facePatchTopo
+            (    
+                markFacesOnProblemCells
+                (
+                    globalToPatch
+                )
+            );
+            forAll(facePatchTopo, faceI)
+            {
+                if (facePatchTopo[faceI] != -1)
+                {
+                    problemTopo.insert(faceI);
+                }
+            }
+            Pout<< "Dumping " << problemTopo.size()
+                << " problem faces to " << problemTopo.objectPath() << endl;
+            problemTopo.write();
+        }
+
         Info<< "Introducing baffles to delete problem cells." << nl << endl;
 
         if (debug)
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index b8ee73b8ce..e9eb2cb445 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -1277,22 +1277,13 @@ Foam::autoPtr<Foam::mapDistributePolyMesh>
 
     if (Pstream::nProcs() > 1)
     {
-        labelList distribution(decomposer.decompose(mesh_.cellCentres()));
-        // Get distribution such that baffle faces stay internal to the
-        // processor.
-        //labelList distribution(decomposePreserveBaffles(decomposer));
-
-        if (debug)
-        {
-            Pout<< "Wanted distribution:"
-                << distributor.countCells(distribution)
-                << endl;
-        }
-        // Do actual sending/receiving of mesh
-        distMap = distributor.distribute(distribution);
-
-        // Update numbering of meshRefiner
-        distribute(distMap);
+        distMap = balance
+        (
+            false,  //keepZoneFaces
+            false,  //keepBaffles
+            decomposer,
+            distributor
+        );
 
         Info<< "Balanced mesh in = "
             << mesh_.time().cpuTimeIncrement() << " s" << endl;
@@ -1302,7 +1293,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh>
 
         if (debug)
         {
-            Pout<< "Writing " << msg
+            Pout<< "Writing balanced " << msg
                 << " mesh to time " << mesh_.time().timeName() << endl;
             write
             (
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index ce09e597e2..1654e3cd7f 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -377,15 +377,34 @@ Foam::refinementSurfaces::refinementSurfaces
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Get indices of named surfaces (surfaces with cellZoneName)
+// Get indices of unnamed surfaces (surfaces without faceZoneName)
+Foam::labelList Foam::refinementSurfaces::getUnnamedSurfaces() const
+{
+   labelList anonymousSurfaces(faceZoneNames_.size());
+
+    label i = 0;
+    forAll(faceZoneNames_, surfI)
+    {
+        if (faceZoneNames_[surfI].size() == 0)
+        {
+            anonymousSurfaces[i++] = surfI;
+        }
+    }
+    anonymousSurfaces.setSize(i);
+
+    return anonymousSurfaces;
+}
+
+
+// Get indices of named surfaces (surfaces with faceZoneName)
 Foam::labelList Foam::refinementSurfaces::getNamedSurfaces() const
 {
-   labelList namedSurfaces(cellZoneNames_.size());
+   labelList namedSurfaces(faceZoneNames_.size());
 
     label namedI = 0;
-    forAll(cellZoneNames_, surfI)
+    forAll(faceZoneNames_, surfI)
     {
-        if (cellZoneNames_[surfI].size() > 0)
+        if (faceZoneNames_[surfI].size() > 0)
         {
             namedSurfaces[namedI++] = surfI;
         }
@@ -846,6 +865,69 @@ void Foam::refinementSurfaces::findNearest
 }
 
 
+void Foam::refinementSurfaces::findNearestRegion
+(
+    const labelList& surfacesToTest,
+    const pointField& samples,
+    const  scalarField& nearestDistSqr,
+    labelList& hitSurface,
+    labelList& hitRegion
+) const
+{
+    labelList geometries(IndirectList<label>(surfaces_, surfacesToTest));
+
+    // Do the tests. Note that findNearest returns index in geometries.
+    List<pointIndexHit> hitInfo;
+    searchableSurfacesQueries::findNearest
+    (
+        allGeometry_,
+        geometries,
+        samples,
+        nearestDistSqr,
+        hitSurface,
+        hitInfo
+    );
+
+    // Rework the hitSurface to be surface (i.e. index into surfaces_)
+    forAll(hitSurface, pointI)
+    {
+        if (hitSurface[pointI] != -1)
+        {
+            hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
+        }
+    }
+
+    // Collect the region
+    hitRegion.setSize(hitSurface.size());
+    hitRegion = -1;
+
+    forAll(surfacesToTest, i)
+    {
+        label surfI = surfacesToTest[i];
+
+        // Collect hits for surfI
+        const labelList localIndices(findIndices(hitSurface, surfI));
+
+        List<pointIndexHit> localHits
+        (
+            IndirectList<pointIndexHit> 
+            (
+                hitInfo,
+                localIndices
+            )
+        );
+
+        labelList localRegion;
+        allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
+
+        forAll(localIndices, i)
+        {
+            hitRegion[localIndices[i]] = localRegion[i];
+        }
+    }
+}
+
+
 //// Find intersection with max of edge. Return -1 or the surface
 //// with the highest maxLevel above currentLevel
 //Foam::label Foam::refinementSurfaces::findHighestIntersection
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index edb613d966..ee5dd219a5 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -154,7 +154,10 @@ public:
                 return cellZoneNames_;
             }
 
-            //- Get indices of named surfaces (surfaces with cellZoneName)
+            //- Get indices of named surfaces (surfaces with faceZoneName)
+            labelList getUnnamedSurfaces() const;
+
+            //- Get indices of named surfaces (surfaces with faceZoneName)
             labelList getNamedSurfaces() const;
 
             //- Get indices of closed named surfaces
@@ -249,6 +252,8 @@ public:
 
             //- Find intersection nearest to the endpoints. surface1,2 are
             //  not indices into surfacesToTest but refinement surface indices.
+            //  Returns surface, region on surface (so not global surface)
+            //  and position on surface.
             void findNearestIntersection
             (
                 const labelList& surfacesToTest,
@@ -282,6 +287,17 @@ public:
                 List<pointIndexHit>&
             ) const;
 
+            //- Find nearest point on surfaces. Return surface and region on
+            //  surface (so not global surface)
+            void findNearestRegion
+            (
+                const labelList& surfacesToTest,
+                const pointField& samples,
+                const  scalarField& nearestDistSqr,
+                labelList& hitSurface,
+                labelList& hitRegion
+            ) const;
+
             //- Detect if a point is 'inside' (closed) surfaces.
             //  Returns -1 if not, returns first surface it is.
             void findInside
diff --git a/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C b/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
index c5927a6128..907eb097f1 100644
--- a/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
+++ b/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
@@ -301,7 +301,7 @@ void Foam::layerAdditionRemoval::removeCellLayer
 
         // Is any of the faces a boundary face?  If so, grab the patch
         // A boundary-to-boundary collapse is checked for in validCollapse()
-        // and cannnot happen here.  
+        // and cannot happen here.  
 
         if (!mesh.isInternalFace(mf[faceI]))
         {
diff --git a/src/meshTools/searchableSurface/searchableBox.C b/src/meshTools/searchableSurface/searchableBox.C
index b1776b3c9e..1858905d49 100644
--- a/src/meshTools/searchableSurface/searchableBox.C
+++ b/src/meshTools/searchableSurface/searchableBox.C
@@ -436,7 +436,6 @@ void Foam::searchableBox::findLineAll
     // Work array
     DynamicList<pointIndexHit, 1, 1> hits;
 
-//XXX
     // Tolerances:
     // To find all intersections we add a small vector to the last intersection
     // This is chosen such that
diff --git a/tutorials/pimpleFoam/t-junction/0/U b/tutorials/pimpleFoam/t-junction/0/U
new file mode 100644
index 0000000000..8851e66583
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/U
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet      
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+
+    outlet1    
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
+    }
+
+    outlet2   
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
+    }
+
+    defaultFaces      
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/0/epsilon b/tutorials/pimpleFoam/t-junction/0/epsilon
new file mode 100644
index 0000000000..8e5b3e553d
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/epsilon
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    inlet
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.01;         // 1cm - half channel height
+        value           uniform 1;
+    }
+
+    outlet1
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+    }
+
+    outlet2
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+    }
+
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/0/k b/tutorials/pimpleFoam/t-junction/0/k
new file mode 100644
index 0000000000..35b2bd3ff9
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/k
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;           // 5% turbulence
+        value           uniform 1;
+    }
+
+    outlet1
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+    }
+
+    outlet2
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+    }
+
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/0/nuTilda b/tutorials/pimpleFoam/t-junction/0/nuTilda
new file mode 100644
index 0000000000..6ea76a6b63
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/nuTilda
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nuTilda;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet1
+    {
+        type            zeroGradient;
+    }
+
+    outlet2
+    {
+        type            zeroGradient;
+    }
+
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/0/p b/tutorials/pimpleFoam/t-junction/0/p
new file mode 100644
index 0000000000..b3d4646058
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/p
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    inlet      
+    {
+        type            totalPressure;
+        p0              uniform 100040;
+        U               U;
+        phi             phi;
+        rho             none;
+        psi             none;
+        gamma           1;
+        value           uniform 100040;
+    }
+
+    outlet1     
+    {
+        type            fixedValue;
+        value           uniform 100010;
+    }
+
+    outlet2    
+    {
+        type            fixedValue;
+        value           uniform 100000;
+    }
+
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/README.txt b/tutorials/pimpleFoam/t-junction/README.txt
new file mode 100644
index 0000000000..0905f8784d
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/README.txt
@@ -0,0 +1,2 @@
+15/8/8 Simple T-junction. Inlet on left, one outlet at bottom, one at top.
+To test multiple outlets.
diff --git a/tutorials/pimpleFoam/t-junction/constant/RASProperties b/tutorials/pimpleFoam/t-junction/constant/RASProperties
new file mode 100644
index 0000000000..e46e0a858b
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/constant/RASProperties
@@ -0,0 +1,200 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel            kEpsilon;
+
+turbulence          on;
+
+printCoeffs         on;
+
+laminarCoeffs
+{
+}
+
+kEpsilonCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphaEps         0.76923;
+}
+
+RNGkEpsilonCoeffs
+{
+    Cmu              0.0845;
+    C1               1.42;
+    C2               1.68;
+    alphak           1.39;
+    alphaEps         1.39;
+    eta0             4.38;
+    beta             0.012;
+}
+
+realizableKECoeffs
+{
+    Cmu              0.09;
+    A0               4.0;
+    C2               1.9;
+    alphak           1;
+    alphaEps         0.833333;
+}
+
+kOmegaSSTCoeffs
+{
+    alphaK1          0.85034;
+    alphaK2          1.0;
+    alphaOmega1      0.5;
+    alphaOmega2      0.85616;
+    gamma1           0.5532;
+    gamma2           0.4403;
+    beta1            0.0750;
+    beta2            0.0828;
+    betaStar         0.09;
+    a1               0.31;
+    c1               10;
+
+    Cmu              0.09;
+}
+
+NonlinearKEShihCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphak           1;
+    alphaEps         0.76932;
+    A1               1.25;
+    A2               1000;
+    Ctau1            -4;
+    Ctau2            13;
+    Ctau3            -2;
+    alphaKsi         0.9;
+}
+
+LienCubicKECoeffs
+{
+    C1               1.44;
+    C2               1.92;
+    alphak           1;
+    alphaEps         0.76923;
+    A1               1.25;
+    A2               1000;
+    Ctau1            -4;
+    Ctau2            13;
+    Ctau3            -2;
+    alphaKsi         0.9;
+}
+
+QZetaCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphaZeta        0.76923;
+    anisotropic     no;
+}
+
+LaunderSharmaKECoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphaEps         0.76923;
+}
+
+LamBremhorstKECoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphaEps         0.76923;
+}
+
+LienCubicKELowReCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphak           1;
+    alphaEps         0.76923;
+    A1               1.25;
+    A2               1000;
+    Ctau1            -4;
+    Ctau2            13;
+    Ctau3            -2;
+    alphaKsi         0.9;
+    Am               0.016;
+    Aepsilon         0.263;
+    Amu              0.00222;
+}
+
+LienLeschzinerLowReCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphak           1;
+    alphaEps         0.76923;
+    Am               0.016;
+    Aepsilon         0.263;
+    Amu              0.00222;
+}
+
+LRRCoeffs
+{
+    Cmu              0.09;
+    Clrr1            1.8;
+    Clrr2            0.6;
+    C1               1.44;
+    C2               1.92;
+    Cs               0.25;
+    Ceps             0.15;
+    alphaEps         0.76923;
+}
+
+LaunderGibsonRSTMCoeffs
+{
+    Cmu              0.09;
+    Clg1             1.8;
+    Clg2             0.6;
+    C1               1.44;
+    C2               1.92;
+    C1Ref            0.5;
+    C2Ref            0.3;
+    Cs               0.25;
+    Ceps             0.15;
+    alphaEps         0.76923;
+    alphaR           1.22;
+}
+
+SpalartAllmarasCoeffs
+{
+    alphaNut         1.5;
+    Cb1              0.1355;
+    Cb2              0.622;
+    Cw2              0.3;
+    Cw3              2;
+    Cv1              7.1;
+    Cv2              5.0;
+}
+
+wallFunctionCoeffs
+{
+    kappa            0.4187;
+    E                9;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/constant/polyMesh/blockMeshDict b/tutorials/pimpleFoam/t-junction/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000..d8f6e3bb6e
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/constant/polyMesh/blockMeshDict
@@ -0,0 +1,109 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//           outlet1
+//             +-+
+//             | |
+//             | |
+//             | |
+//             | |
+// +-----------+ |
+// |inlet        |
+// +-----------+ |
+//             | |
+//             | |
+//             | |
+//             | |
+//             +-+
+//           outlet2
+
+convertToMeters 1;
+
+vertices        
+(
+    (0.0  -0.01 0)   //0
+    (0.2  -0.01 0)
+    (0.2   0.01 0)   //2
+    (0.0   0.01 0)
+
+    (0.22 -0.01 0)  //4
+    (0.22  0.01 0)
+
+    (0.2  -0.21 0)  //6
+    (0.22 -0.21 0)
+
+    (0.2   0.21 0)  //8
+    (0.22  0.21 0)
+
+    // Z
+    (0.0  -0.01 0.02)   //0
+    (0.2  -0.01 0.02)
+    (0.2   0.01 0.02)   //2
+    (0.0   0.01 0.02)
+
+    (0.22 -0.01 0.02)  //4
+    (0.22  0.01 0.02)
+
+    (0.2  -0.21 0.02)  //6
+    (0.22 -0.21 0.02)
+
+    (0.2   0.21 0.02)  //8
+    (0.22  0.21 0.02)
+
+);
+
+blocks          
+(
+    // inlet block
+    hex (0 1 2 3  10 11 12 13) (50 5 5) simpleGrading (1 1 1)
+
+    // central block
+    hex (1 4 5 2  11 14 15 12) (5 5 5) simpleGrading (1 1 1)
+
+    // bottom block
+    hex (6 7 4 1  16 17 14 11) (5 50 5) simpleGrading (1 1 1)
+
+    // top block
+    hex (2 5 9 8  12 15 19 18) (5 50 5) simpleGrading (1 1 1)
+);
+
+edges           
+(
+);
+
+patches         
+(
+    patch inlet
+    (
+        (0 10 13 3)
+    )
+
+    patch outlet1
+    (
+        (6 7 17 16)
+    )
+
+    patch outlet2
+    (
+        (8 18 19 9)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/constant/polyMesh/boundary b/tutorials/pimpleFoam/t-junction/constant/polyMesh/boundary
new file mode 100644
index 0000000000..4ea7f95a3e
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/constant/polyMesh/boundary
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+4
+(
+    inlet
+    {
+        type            patch;
+        nFaces          25;
+        startFace       10050;
+    }
+    outlet1
+    {
+        type            patch;
+        nFaces          25;
+        startFace       10075;
+    }
+    outlet2
+    {
+        type            patch;
+        nFaces          25;
+        startFace       10100;
+    }
+    defaultFaces
+    {
+        type            wall;
+        nFaces          3075;
+        startFace       10125;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/constant/transportProperties b/tutorials/pimpleFoam/t-junction/constant/transportProperties
new file mode 100644
index 0000000000..fd5fc94660
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/constant/transportProperties
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              nu [0 2 -1 0 0 0 0] 1e-05;
+
+CrossPowerLawCoeffs
+{
+    nu0             nu0 [0 2 -1 0 0 0 0] 1e-06;
+    nuInf           nuInf [0 2 -1 0 0 0 0] 1e-06;
+    m               m [0 0 1 0 0 0 0] 1;
+    n               n [0 0 0 0 0 0 0] 1;
+}
+
+BirdCarreauCoeffs
+{
+    nu0             nu0 [0 2 -1 0 0 0 0] 1e-06;
+    nuInf           nuInf [0 2 -1 0 0 0 0] 1e-06;
+    k               k [0 0 1 0 0 0 0] 0;
+    n               n [0 0 0 0 0 0 0] 1;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/system/controlDict b/tutorials/pimpleFoam/t-junction/system/controlDict
new file mode 100644
index 0000000000..b5f0d20f7e
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/system/controlDict
@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application turbFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+maxCo           5;
+
+functions
+(
+    probes
+    {
+        // Type of functionObject
+        type probes;
+
+        // Where to load it from (if not already in solver)
+        functionObjectLibs ("libsampling.so");
+
+        // Name of the directory for the probe data
+        name probes;
+
+        // Locations to be probed. runTime modifiable!
+        probeLocations
+        (
+            (1E-6 0 0.01)           // at inlet
+            (0.21 -0.20999 0.01)    // at outlet1
+            (0.21  0.20999 0.01)    // at outlet2
+            (0.21 0 0.01)           // at central block
+        );
+
+        // Fields to be probed. runTime modifiable!
+        fields
+        (
+            p
+            U
+        );
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/system/fvSchemes b/tutorials/pimpleFoam/t-junction/system/fvSchemes
new file mode 100644
index 0000000000..bbe7fab0d2
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/system/fvSchemes
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+    grad(U)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss limitedLinearV 1;
+    div(phi,k)      Gauss limitedLinear 1;
+    div(phi,epsilon) Gauss limitedLinear 1;
+    div(phi,R)      Gauss limitedLinear 1;
+    div(R)          Gauss linear;
+    div(phi,nuTilda) Gauss limitedLinear 1;
+    div((nuEff*dev(grad(U).T()))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(nuEff,U) Gauss linear corrected;
+    laplacian((1|A(U)),p) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DREff,R) Gauss linear corrected;
+    laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(U)  linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/system/fvSolution b/tutorials/pimpleFoam/t-junction/system/fvSolution
new file mode 100644
index 0000000000..995edf623c
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/system/fvSolution
@@ -0,0 +1,88 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p GAMG
+    {
+        tolerance        1e-6;
+        relTol           0.01;
+
+        smoother         GaussSeidel;
+
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+    };
+    pFinal GAMG
+    {
+        tolerance        1e-6;
+        relTol           0.0;
+
+        smoother         GaussSeidel;
+
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+    };
+
+    U PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0.1;
+    };
+    UFinal PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    k PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    epsilon PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+}
+
+PIMPLE
+{
+    nOuterCorrectors 2;
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+relaxationFactors
+{
+    //p               0.3;
+    U               1.0;
+    k               1.0;
+    epsilon         1.0;
+}
+
+
+// ************************************************************************* //